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Introduction 

The  expected  long-term  outcome  of  this  research  is  the  enhancement  of  microwave  breast 
cancer  detection  sensitivity  through  the  use  of  extrinsic  contrast  mechanisms,  namely 
biocompatible  micro-  or  nano-particles  systemically  delivered  to  the  site  of  a  malignancy  via 
passive  or  active  means.  Such  particles  may  impact  the  malignant-to-normal  contrast  of  breast 
tissue  at  microwave  frequencies,  thereby  altering  the  microwave  scattering,  absorption,  and 
thermoacoustic  response  of  breast  tumors  and  enhancing  the  overall  efficacy  of  microwave 
breast  cancer  detection.  The  goal  of  this  project  was  to  conduct  a  series  of  experimental  and 
numerical  investigations  to  elucidate  and  optimize  the  microwave-frequency  response  of 
microparticle  and  nanoparticle  contrast  agents.  The  specific  aims  were  as  follows:  (1) 
Experimentally  characterize  the  effective  electromagnetic  properties  of  a  variety  of 
micro/nanoparticle  dispersions  at  microwave  frequencies;  (2)  Develop  theoretical/numerical 
tools  to  model  microwave  interactions  with  micro/nanoparticles  and  utilize  these  tools  to  design 
micro/nanoparticles  that  optimize  the  microwave  contrast  of  malignant  breast  tumors;  and  (3) 
Conduct  experimental  investigations  of  the  microwave  scattering  and  absorption  properties  of 
the  micro/nanoparticle  contrast  agents  in  tissue-mimicking  phantoms.  These  three  aims  were 
addressed  as  three  major  tasks  described  below.  The  discussion  below  refers  to  results  in 
published  papers  and  abstracts  included  in  the  Appendices.  Any  figure  that  is  available  in  an 
appendix  is  not  duplicated  in  the  body  of  this  report. 

Body 

Task  1.  Experimentally  characterize  the  effective  electromagnetic  properties  of  a  variety  of 
micro/nanoparticle  dispersions  at  microwave  frequencies. 

Task  1.1.  Our  studies  focused  on  the  subset  of  candidate  micro-  and  nano-particles  that  are 
the  most  promising  for  altering  the  scattering  and  absorption  of  microwave  energy  by  a 
malignant  lesion.  Namely,  we  have  investigated  air-filled  microbubbles,  metallic  nanoparticles, 
and  carbon  nanotubes.  The  goal  of  the  experimental  studies  of  Task  1  was  to  characterize  the 
extent  to  which  the  dielectric  properties  of  particles-in-liquid  mixtures  differed  from  that  of  the 
background  liquid  alone,  without  the  particles.  Since  there  was  no  biocompatibility  requirement 
for  these  laboratory  studies,  we  had  some  flexibility  in  choosing  the  tissue-mimicking 
background  liquid.  The  liquid  dispersions  that  we  acquired  or  created  involved  either  ethylene 
glycol  (EG)  or  water  as  the  background  liquid.  These  were  chosen  as  inexpensive,  readily 
available  reference  liquids  with  well  characterized  dielectric  properties  at  microwave 
frequencies. 

•  Microbubbles:  We  acquired  air-filled,  glass  spheres  from  3M  (iM30K,  average  diameter  of 
18  microns)  and  dispersed  them  in  EG  to  create  mixtures  of  various  concentrations. 

•  Metallic  nanoparticles:  We  acquired  monodispersions  of  metallic  nanoparticles  from  a) 
Meliorum  (copper  nanoparticles  dispersed  in  EG,  with  an  average  diameter  of  384  nm  and  an 
overall  volume  fraction  of  30%)  and  b)  the  Messersmith  laboratory  at  Northwestern 
University  (spherical  and  rod-shaped  gold  nanoparticles  synthesized  using  an  aqueous  seeded 
growth  method  and  coated  with  CTAB  surfactant,  with  an  average  diameter  of  20  nm  and  in 
concentrations  ranging  from  0. 1  to  20  mM). 
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•  Carbon  nanotubes:  We  acquired  monodispersions  of  single-walled  carbon  nanotubes 
(SWCNTs)  dispersed  in  aqueous  solutions  of  1%  Pluronic  (F127)  surfactant  from  Prof. 
Balaji  Sitharaman  in  the  Department  of  Biomedical  Engineering  at  Stony  Brook  University. 
We  also  acquired  SWCNTs  synthesized  by  commercial  HiPCO  (Unidym)  and  functionalized 
with  a  biocompatible  and  stable  surfactant  (polyethylene-glycol(PEG)-ylated  phospholipids 
(PL)). 

Task  1.2.  Prior  to  characterizing  the  dielectric  properties  of  the  micro/nanoparticle 
dispersions  and  the  corresponding  particle-free  liquids,  we  conducted  microwave-frequency 
measurements  of  the  dielectric  properties  of  well-known  homogeneous  reference  liquids  (without 
particles)  as  a  validation  step.  We  performed  the  measurements  using  two  complementary 
techniques:  a  wideband  (0.5-20  GHz)  open-ended  coaxial  probe  technique  that  we  previously 
perfected  and  used  in  a  large-scale  study  of  breast  tissue  [1-4],  and  a  narrowband  cavity 
perturbation  technique  [5]  conducted  at  2.2,  5.05,  and  8  GHz.  Figure  1  shows  a  comparison  of 
the  measured  dielectric  properties  of  water,  methanol,  ethanol,  and  EG  obtained  using  these  two 
characterization  techniques.  The  excellent  agreement  observed  between  the  narrowband  and 
wideband  results,  as  well  as  the  agreement  with  reported  values  in  the  literature,  validates  our 
approach  and  illustrates  the  consistency  check  offered  by  the  two  complementary  techniques. 


Figure  1.  Dielectric  properties  (relative  permittivity,  left,  and  effective  conductivity,  right)  of  four 
reference  liquids  measured  using  a  wideband  technique  (solid  lines)  and  a  narrowband  technique  (data 
markers). 

Task  1.3.  We  used  the  wideband  and  narrowband  microwave-frequency  measurement 
techniques  described  above  to  characterize  the  effective  dielectric  properties  of  the  three  classes 
of  micro/nanoparticle  dispersions.  We  also  analyzed  the  measured  data  to  assess  changes  in  the 
electromagnetic  properties  of  the  dispersion  liquid  due  to  the  presence  of  the 
micro/nanoparticles.  A  summary  of  our  findings  from  this  task  is  provided  below;  further  details 
about  the  experimental  studies  are  provided  in  Appendices  A  and  E. 

Task  1.3a:  Microbubbles.  The  reprint  provided  in  Appendix  A  (Figure  2  therein)  reports  the 
relative  permittivity  and  effective  conductivity  of  pure  EG  (Solution  1)  and  four  different 
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microbubble-EG  dispersions  (Solutions  2-5).  The  microbubble  concentrations,  ranging  from 
20%  (Solution  2)  up  to  40%  by  weight  (Solution  5),  are  within  the  typical  109  microbubbles/mL 
dosage  administered  in  human  ultrasound  procedures  [6],  We  found  that  the  dielectric  properties 
decrease  with  increasing  microbubble  concentration,  as  expected.  In  particular,  the  relative 
permittivity  and  effective  conductivity  at  3  GHz  of  the  highest-concentration  dispersion 
decreased  by  factors  of  approximately  two  and  three,  respectively.  Thus,  the  measured  change  in 
dielectric  properties  is  significant  for  clinically  relevant  concentrations. 


Task  1.3b:  Metallic  nanoparticles.  Figure  2  below  shows  the  relative  permittivity  and 
effective  conductivity  of  the  commercial  dispersion  consisting  of  a  30%  volume  fraction  of  the 
copper  nanoparticles  in  EG.  We  found  reasonable  agreement  between  the  wideband  (solid  and 
dotted  lines)  and  narrowband  measurements  (symbols).  For  both  the  probe  and  cavity 
measurements  above  2  GHz,  the  relative  permittivity  (dielectric  constant)  of  the  dispersion  is 
approximately  20%  higher  than  that  of  pure  EG.  However,  we  observed  no  significant  change  in 
the  conductivity  due  to  the  presence  of  the  nanoparticles.  We  also  characterized  the  dielectric 
properties  of  the  Northwestern  University  samples  of  gold  nanoparticles  in  an  aqueous  solution 
and  found  no  statistically  significant  change  relative  to  the  particle-free  background  medium. 
These  negative  findings  motivated  our  investigation  of  SWCNTs  as  high-aspect-ratio  alternative 
particles  to  compact  metallic  nanoparticles,  as  discussed  further  below. 


Figure  2.  Dielectric  properties  of  a  commercially  prepared  dispersion  of  384-nm  copper  particles  in  EG 
(30%  volume  fraction)  compared  with  those  of  pure  EG  (the  background  medium).  The  measurements 
were  conducted  using  both  a  wideband  open-ended  coaxial  probe  technique  ("Probe")  and  a  narrowband 
cavity  perturbation  technique  ("Cavity"). 

Task  1.3c:  SWCNTs.  We  characterized  the  dielectric  properties  of  the  Stony  Brook  CNTs 
dispersed  in  Pluronic  and  water.  First,  we  compared  the  properties  of  monodispersions  consisting 
of  two  concentrations  (0.5  and  1.0  mg/mL)  of  0.5 -micron-long  single-walled  CNTs  with  the 
properties  of  the  pure  aqueous  background.  We  found  that  the  relative  permittivity  and  effective 
conductivity  of  the  CNT  mixtures  are  significantly  higher  than  the  background  medium.  Second, 
we  compared  the  properties  of  monodispersions  consisting  of  two  concentrations  (0.5  and  1.0 
mg/mL)  of  buckyballs  with  the  properties  of  the  pure  aqueous  background.  In  contrast  to  the 
large  change  observed  with  the  CNTs,  we  found  that  the  buckyballs  did  not  change  the  effective 
properties  of  the  medium.  This  finding  is  consistent  with  our  measurement  results  for  the 
metallic  nanoparticles,  which  are  similar  to  the  buckyballs  in  that  they  are  compact  in  shape  with 
small  aspect  ratios.  Thus  our  results  indicate  that  for  relatively  small  volume  fractions  of 
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nanoparticles,  large  aspect  ratios  are  required  in  the  metallic  or  semi-metal  particles  in  order  to 
significantly  alter  the  macroscopic  effective  dielectric  properties  of  the  medium. 

We  also  characterized  the  dielectric  properties  of  dispersions  created  from  the  biocompatible 
(PL-PEG)  functionalized  SWCNTs  in  water  (see  Appendix  E).  Key  results  of  this  study  are 
summarized  in  Figure  3.  For  a  concentration  of  2  mg/mL  of  SWCNTs,  we  observed  at  1  GHz  an 
increase  in  relative  permittivity  of  approximately  10%  and  an  increase  in  effective  conductivity 
of  approximately  90%  relative  to  the  properties  of  PL-PEG- water  only  (no  SWCNTs). 


Figure  3.  Frequency-dependent  relative  permittivity  and  effective  conductivity  of  a  mixture  of  water  and 
phospholipid  (PL)-polyethylene-glycol  (PEG)  with  and  without  SWCNTs  at  a  concentration  of  2  mg/mL. 
The  measurements  were  conducted  using  a  precision  open-ended  coaxial  probe  technique. 

Task  2.  Develop  theoretical/numerical  tools  to  model  microwave  interactions  with 
micro/nanoparticles  and  utilize  these  tools  to  design  micro/nanoparticles  that  optimize  the 
microwave  contrast  of  malignant  breast  tumors. 

Task  2.1.  We  adapted  two  types  of  2D  classical  electromagnetic  simulation  tools  and 
conducted  a  comprehensive  theoretical  investigation  of  the  effects  of  moderate  volume  fractions 
of  a  variety  of  micro/nanoparticles  on  the  effective  electromagnetic  properties  of  binary  mixtures 
with  a  background  medium  representative  of  malignant  breast  tissue.  Since  2D  simulations 
require  only  a  fraction  of  the  computational  resources  of  3D  simulations,  we  used  2D  methods  to 
perform  a  large  number  (1000's)  of  simulations  with  systematically  varying  particle  sizes  and 
randomized  particle  distributions.  One  of  the  2D  simulation  tools  was  based  on  a  full-wave 
technique  for  solving  Maxwell's  equations  of  electrodynamics,  namely  the  finite-difference  time- 
domain  (FDTD)  method,  while  the  other  tool  was  based  on  a  quasi-static  finite-element  method 
(QS-FEM).  The  computational  burden  of  the  latter  technique  is  lower,  but  its  accuracy  is  limited 
in  the  case  of  lossy  particles  as  it  does  not  account  for  physical  size  effects  such  as  skin  depth. 

We  focused  on  lossless  air- filled  microbubbles  and  conducting  metallic  spheres  as  the 
candidate  contrast  agent  particles  in  this  study.  We  modeled  particles  that  ranged  in  size  from 
500  nm  to  200  microns.  We  expect  clinically  relevant  microbubble  contrast  agents  to  be  1-7 
microns  in  diameter,  and  metallic  spheres  to  range  in  size  from  200  to  250  nm.  While  the 
smallest  metallic  particles  considered  in  this  numerical  study  exceed  the  predicted  size  range  for 
practical  applications,  we  limited  the  scope  of  this  2D  classical  electromagnetics  study  to  a  size 
regime  that  poses  the  least  uncertainty  about  the  physical  properties  of  the  metallic  particles. 
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Furthermore,  we  chose  the  largest  particle  size  to  be  approximately  an  order  of  magnitude  larger 
than  the  skin  depth  at  microwave  frequencies  in  order  to  investigate  the  effects  of  inclusion  size 
on  the  effective  dielectric  properties  of  the  mixture.  For  both  the  QS-FEM  and  FDTD 
simulations,  we  considered  5%,  10%,  and  20%  volume  fractions,  which  we  estimate  to  be 
reasonable  for  clinical  applications.  For  each  volume  fraction,  we  simulated  100  spatial 
arrangements  of  particles  and  computed  ensemble  averages  of  the  effective  dielectric  properties. 
We  summarize  our  findings  below;  further  details  about  the  2D  numerical  studies  of 
microbubbles  and  metallic  nanoparticle  dispersions  are  provided  in  Appendix  D. 

We  found  excellent  agreement  between  the  QS-FEM-computed  and  FDTD-computed 
effective  dielectric  properties  of  mixtures  containing  microbubbles,  as  expected  for  the  case  of 
lossless  particles  where  the  simplifying  assumptions  of  the  QS-FEM  are  valid.  We  also  found 
excellent  agreement  between  the  numerical  results  obtained  here  for  the  microbubbles  and  the 
experimental  results  obtained  under  Task  1.3  that  suggested  that  small-to-moderate  volume 
fractions  of  microbubbles  can  substantially  lower  the  effective  dielectric  properties  of  the 
mixture. 

For  the  metallic  particle  simulations,  we  found  great  discrepancies  between  the  QS-FEM  and 
FDTD  values,  and  considered  this  to  be  evidence  that  the  quasi-static  assumption  was  breaking 
down  for  the  lossy  particles.  Therefore  we  relied  upon  the  FDTD  predictions  for  analyzing  the 
effects  of  the  metallic  particles.  We  investigated  the  effect  of  particle  spatial  distribution  (e.g. 
contacting  versus  non-contacting)  and  particle  size  on  the  effective  dielectric  properties  of  the 
mixtures.  We  found  that  contacting  lossy  particles  lead  to  a  larger  increase  in  the  effective 
properties  than  non-contacting  inclusions.  We  observed  an  interesting  size-dependent  trend  in 
the  effective  properties  as  the  inclusion  diameter  increased  from  one-tenth  to  ten  times  the  skin 
depth  at  the  microwave  frequency  assumed  in  the  simulations  (5  GHz).  Specifically,  we  found 
that  when  the  inclusion  size  is  approximately  one  to  ten  times  the  skin  depth,  the  relative 
permittivity  of  the  mixture  decreases  rapidly  while  the  effective  conductivity  reaches  a 
maximum. 

We  also  developed  efficient  3-D  algorithms  to  calculate  the  microwave  electromagnetic 
properties  of  binary  mixtures  of  micro/nanoparticles.  Our  algorithm  development  was  based  on 
the  generalized  Mie  (GMM)  solution  [7],  which  is  capable  of  solving  Maxwell’s  equations  for 
the  scattering  problem  of  a  particle  ensemble  composed  of  a  large  number  of  spherical  particles. 
Our  original  plan  was  to  benchmark  the  GMM  simulation  results  against  a  3-D  FDTD  model  for 
a  small  number  (N  <  10)  of  particles  and  then  use  it  to  extend  the  2-D  theoretical  investigation 
described  above  to  3-D.  However,  once  our  experimental  results  revealed  a  negligible  dielectric- 
properties  enhancement  associated  with  compact  metallic  nanoparticles,  we  shifted  our  focus 
carbon  nanotube  dispersions. 

Task  2.2.  The  key  first  step  towards  developing  a  first-principles  simulation  tool  for 
modeling  microwave  interactions  with  metallic  nanoparticles  was  to  create  a  hybrid  solver  that 
couples  a  computational  electromagnetics  model  with  a  diffusive  carrier  transport  model. 
Accordingly,  we  developed  a  multi-physics  simulation  tool  that  combines  FDTD  for  solving  the 
time-varying  Maxwell’s  equations  with  an  ensemble  Monte  Carlo  (EMC)  approach  for 
stochastically  solving  the  Boltzman  transport  equation.  We  characterized  the  EMC-FDTD  solver 
in  terms  of  grid  spacing  and  carrier  ensemble  size  (see  Appendix  F).  We  rigorously  validated  this 
state-of-the-art  simulation  tool,  as  described  in  Appendix  B,  by  simulating  the  terahertz- 
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frequency  skin  effect  in  doped  silicon.  We  obtained  excellent  agreement  between  our  EMC- 
FDTD-computed  effective  conductivity  of  silicon  at  low  doping  densities  and  the  measured 
effective  conductivity  reported  by  Jeon  and  Grischkowsky  [8],  The  promising  nature  of  this 
multi-physics  simulation  tool  led  to  more  extensive  development  of  the  EMC-FDTD  solver 
pursued  under  a  separate  project  for  characterizing  the  THz-frequency  conductivity  of  silicon  at 
higher  doping  densities.  The  second  step  of  extending  this  work  to  metal  nanoparticles  was  no 
longer  needed  once  our  focus  shifted  to  CNTs,  as  microwave  interactions  with  CNTs  can  be 
analyzed  rigorously  using  classical/semi-classical  electromagnetic  models. 

Task  2.3.  As  noted  above,  the  microwave  response  of  compact  metallic  nanoparticles  is 
simply  too  weak  to  be  of  practical  value  for  enhancing  the  microwave  scattering  and  absorption 
cross-sections  of  malignant  tumors.  The  fortuitous  experimental  discovery  of  the  strong 
microwave  response  of  CNTs  changed  the  nature  of  this  task  from  one  of  designing  the  optimum 
size  and  shape  of  compact  metallic  nanoparticles  via  simulation  to  one  of  identifying  the  range  of 
CNT  lengths  that  optimize  the  potential  microwave  contrast  of  malignant  tumors.  We  found  that 
the  significant  microwave  contrast  enhancement  provided  by  large-aspect-ratio  SWCNTs  does 
not  depend  greatly  on  the  specific  lengths  of  the  particles,  which  is  affected  by  sonication 
duration.  For  example,  we  compared  the  properties  of  monodispersions  consisting  of  either  1.0- 
micron-long  CNTs  or  0.5 -micron-long  CNTs,  both  at  the  same  concentration  of  1.0  mg/mL,  and 
found  no  significant  difference  in  the  effective  dielectric  properties. 

Task  3.  Conduct  experimental  investigations  of  the  microwave  scattering  and  absorption 
properties  of  the  micro/nanoparticle  contrast  agents  in  tissue-mimicking  phantoms. 

Task  3.1.  We  constructed  tissue-mimicking  (TM)  materials  from  oil-in-gelatin  dispersions 
following  the  procedure  outlined  in  [9].  The  dielectric  properties  of  these  materials  can  be 
customized  to  mimic  the  properties  of  a  variety  of  human  soft  tissues  by  controlling  the 
concentrations  of  the  constituents.  The  results  reported  in  [4]  show  that  a  10%-oil  mixture 
adequately  replicates  the  microwave  properties  of  malignant  breast  tissue  over  our  frequency 
range  of  interest.  Furthermore,  these  materials  are  relatively  inexpensive  to  fabricate  and  posess 
long-term  stability.  Therefore,  we  constructed  10%-oil  TM  materials  with  various 
concentrations  of  SWCNTs  as  “tumor”  phantoms.  We  mixed  1,  2,  and  3  mg/mL  of  SWCNTs 
into  a  1%-by- weight  mixture  of  Pluronic  (FI 27)  and  deionized  water.  Then  we  substituted  these 
mixtures  for  the  pure  deionized  water  in  the  TM  recipe  described  in  [9],  The  resulting  percent- 
by-weight  concentrations  of  SWCNTs  in  the  TM  samples  were  0.07%,  0.15%,  and  0.22%, 
respectively.  For  reference,  we  also  constructed  a  pure  TM  mixture  with  the  same  amount  of 
Pluronic  as  the  other  samples,  but  with  no  SWCNTs.  The  dielectric  properties  of  each  of  the  TM 
samples  were  characterized  using  the  open-ended  coaxial  probe  technique.  At  3  GHz,  we  found 
that  SWCNT  concentrations  as  small  as  0.22%  by  weight  increased  the  relative  permittivity  of 
the  TM  material  by  37%  and  the  effective  conductivity  by  81%.  Further  details  are  provided  in 
Appendix  C  and  G. 

Task  3.2.  The  SWCNTs  used  in  the  investigations  conducted  under  this  task  were 
synthesized  by  a  commercial  vendor  using  a  chemical  vapor  deposition  technique  and  acid 
treated  for  purification.  However,  we  briefly  describe  here  the  state-of-the-art  synthesis 
protocols  we  developed  for  two  other  types  of  nanoparticles:  isotropic  (spherical)  gold  and  silver 
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nanoparticles  and  anisotropic  gold  nanorods.  These  protocols  include  modifying  the  nanoparticle 
surface  properties  to  stably  suspend  them  in  saline  solutions  representing  realistic  physiological 
conditions. 

Task  3.2a:  Spherical  nanoparticles.  We  exploited  the  strong  binding  properties  of  DOPA 
peptides  to  anchor  polymers  onto  material  surfaces  to  form  spherical  gold  and  silver 
nanoparticles.  We  used  self-assembling  polymers  to  exert  control  of  particle  shape,  size  and 
properties.  First,  we  synthesized  a  DOPA4  peptide  in  its  solid  phase,  and  linked  it  to  a  PEG2000 
chain.  The  DOPA4  sequence  was  designed  as  a  multifunctional  molecule.  The  first  function  was 
to  mediate  formation  of  self-assembled  nanoaggregates  through  hydrophobic  DOPA  aggregation, 
in  order  to  yield  DOPA-core/PEG-shell  micelles  that  act  as  templates  for  mineralizing  metallic 
nanoparticles.  10  nm  structures  were  identified  in  both  cryo-TEM  and  dynamic  light  scattering 
(DLS)  techniques  -  strong  evidence  for  micellular  formation.  The  second  function  of  the 
repeating  DOPA  peptide  was  the  reduction  of  gold  and  silver  ions  to  form  metal  nanoparticles. 
Control  over  the  size  of  the  gold  nanoparticles  was  achieved  through  PEG-DOPA4  and  metallic 
salt  concentration.  We  have  characterized  the  synthesized  nanoparticles  with  transmission 
electron  microscopy  (TEM)  and  dynamic  light  scattering.  Figure  4  shows  TEM  micrographs  of 
10  nm  gold  (525  nm  plasmon)  and  17  nm  silver  (415  nm  plasmon)  nanoparticles. 


Figure  4.  TEM 
micrographs  of  (a)  10 
nm  gold  and  (b)  17  nm 
silver  nanoparticles 
formed  with  PEG- 
DOPA4  polymers.  Scale 
bar  =  lOOnm. 


Further  experimental  characterization  of  the  synthesized  nanoparticles  demonstrated  that 
these  polymers  not  only  induced  metallic  nanoparticle  formation,  but  also  stably  anchored  PEG 
to  the  surface  in  saline  conditions,  the  third  vital  function  of  the  DOPA  peptide.  For  example, 
after  four  centrifuge  and  washing  steps,  x-ray  photoelectron  spectroscopy  (XPS)  was  performed 
on  gold  nanoparticles  formed  with  the  PEG-DOPA4  polymer  and  signal  from  both  the  metal  and 
the  polymer  were  detected,  implying  stable  linkage  of  polymer  on  the  nanoparticle  surface. 
Further,  PEG-DOPA4-synthesized  gold  nanoparticles  were  tested  in  saline  to  determine  their 
stability  against  aggregation  in  physiological  conditions,  an  important  requirement  for  clinical 
use  as  well  as  further  experimental  measurements  with  phantom  materials.  No  change  in  the 
dynamic  light  scattering  properties  of  the  nanoparticle  dispersions  was  detected  even  up  to  500 
mM  salt  content,  implying  successful  particle  dispersion,  while  the  conventional  citrate- 
stabilized  gold  NPs  aggregated  at  concentrations  as  low  as  10  mM. 

Task  3.2b:  Nanorods.  Using  an  alternative  synthesis  route,  we  fabricated  gold  nanorods  with 
controlled  aspect  ratios.  In  this  case,  we  used  a  two-step  synthesis/coating  strategy.  First,  we 
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used  the  surfactant  CTAB  as  a  reducing  agent  and  template  to  form  CTAB -coated  nanorods  and 
tuned  the  aspect  ratio  of  the  structures  with  control  over  gold  salt  and  seed  concentrations.  In  a 
second  step,  we  replaced  the  CTAB  coating  with  biocompatible  surface-stabilizing  molecules 
with  anti-fouling  PEG5000-thiol.  Briefly,  CTAB-gold  nanorod  suspensions  were  centrifuged, 
and  excess  CTAB-containing  supernatant  solution  was  discarded.  A  PEG-SH  solution  was  added 
to  the  gold  nanorod  pellet,  and  mixed  for  24  hours.  The  solution  was  centrifuged  and  supernatant 
solution  was  discarded;  this  process  was  repeated  twice.  We  characterized  the  synthesized 
nanorods  using  both  TEM  and  optical  extinction  spectroscopy.  A  TEM  micrograph  and  optical 
extinction  spectra  of  synthesized  gold  nanorods  can  be  seen  in  Figure  5.  The  second  step  of 
coating  replacement  was  demonstrated  to  be  an  essential  step  toward  particle  stability.  Surface 
plasma  analysis  indicated  that  PEG-coated  nanorods  remained  more  stable  towards  aggregation 
in  saline  than  either  the  conventional  citrate  or  CTAB  stabilized  nanorods. 
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Figure  5.  (a)  TEM  micrograph  of  gold  nanorods  synthesized  with  the  CTAB-mediated  method  with  an 
aspect  ratio  of  ~6:1.  (b)  Optical  extinction  spectra  of  gold  nanorods  with  various  aspect  ratios.  The  shift 
of  the  optical  resonance  peaks  indicates  that  the  nanorod  aspect  ratio  can  be  tuned  over  a  wide  range. 

Task  3.3.  We  complemented  the  use  of  simple  experimental  phantoms  in  Task  3.4  by  the  use 
of  complex  (anatomically  realistic)  numerical  breast  phantoms  for  the  microwave  scattering 
studies  conducted  here  under  Task  3.3.  The  results  of  the  dielectric  characterization  experiments 
described  in  Task  3.1  were  used  in  microwave  imaging  studies  in  order  to  determine  the  impact 
of  the  SWCNT-induced  dielectric-properties  changes  on  microwave  breast  cancer  detection  [10]. 
Briefly,  we  applied  microwave  inverse  scattering  techniques  to  reconstruct  images  of  numerical 
breast  phantoms  with  a  compact  malignant  lesion.  Two  images  were  reconstructed:  an  image  of  a 
phantom  with  the  endogenous  dielectric  properties  of  the  tumor  and  another  with  elevated 
dielectric  properties  values  due  to  the  carbon  nanotubes.  A  differential  image  was  formed  by 
subtracting  the  two  images  obtained  from  pre-  and  postcontrast  measurements.  Using  this 
differential  imaging  technique,  Shea  et  al.  [10]  reported  detection  of  previously  undetectable 
tumors  in  four  different  classes  of  numerical  phantoms  that  ranged  from  “mostly  fatty”  to 
“extremely  dense”  in  radiographic  density.  These  results  suggest  that  the  changes  in  dielectric 
properties  reported  under  Task  3.1  (and  in  Appendix  C)  are  significant  enough  to  dramatically 
improve  the  sensitivity  of  microwave  imaging.  We  also  created  differential  images  for  the  case 
of  microbubbles  as  the  assumed  contrast  agent  and  obtained  similarly  promising  results  [10]. 
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Task  3.4a:  Microbubbles.  We  characterized  the  thermoacoustic  response  of  simple  “tumors” 
comprised  of  6-mm-diameter  plastic  tubes  fdled  with  microbubble-EG  dispersions.  Appendix  A 
reports  the  details  of  the  experimental  set-up  and  shows  the  thermoacoustic  temporal  waveforms 
measured  for  each  of  six  targets.  The  amplitude  of  the  response  diminishes  and  the  temporal 
width  decreases  with  increasing  concentration  of  microbubbles.  We  conducted  an  analysis  of  the 
effect  of  acoustic  scattering  to  confirm  that  the  observed  reduction  in  the  first  peak  of  the 
thermoacoustic  waveform  is  indeed  due  to  the  reduced  effective  conductivity  (as  shown  in 
Figure  2  of  Appendix  A)  rather  than  due  to  increased  acoustic  attenuation.  With  organic-shell 
microbubbles,  we  expect  the  same  trend  of  reduced  amplitude  but  a  broadening  of  the  temporal 
width. 

Task  3.4b:  Metallic  nanoparticles.  We  conducted  microwave  heating  experiments  on  the 
low-concentration  aqueous  solutions  of  metallic  nanoparticles  synthesized  in  Task  3.2  to 
characterize  the  absorption  and  thermal  properties  of  these  particles.  The  schematic  diagram  of 
the  experimental  setup  is  demonstrated  in  Figure  6.  Briefly,  a  4-GHz  microwave  signal  was 
generated  and  amplified  to  8W.  The  microwave  source  was  matched  to  a  section  of  waveguide  to 
heat  the  nanoparticle  solution  (contained  in  a  plastic  vial)  placed  at  the  waveguide  aperture.  The 
voltage  standing  wave  ratio  (VSWR)  at  the  waveguide  input  was  carefully  calibrated  for  all 
measurements.  A  computer  was  used  to  accurately  control  the  microwave  source  and  the 
heating/cooling  phases  of  the  sample.  A  Luxtron  1652  fiber-optic  temperature  sensor  was 
inserted  to  the  sample  vial  to  record  temperature  in  real  time.  For  each  sample,  the 
heating/cooling  cycles  were  repeated  for  several  times  to  confirm  repeatability.  Figure  7  shows 
the  temperature  trace  of  a  5mM  gold  nanorod  solution  compared  to  deionized  (DI)  water  for 
three  heating/cooling  cycles.  Initially,  we  assumed  that  the  compact  nanoparticles  were 
responsible  for  the  enhanced  heating  response.  However,  we  conducted  additional  experiments 
with  CTAB  surfactant  added  to  the  DI  water  (e.g.,  no  nanoparticles)  and  discovered  that  it  was 
the  CTAB  surfactant  -  not  the  nanoparticles  -  that  was  responsible  for  the  heating  enhancement. 
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Figure  6.  Experimental  setup  for  characterizing  the  microwave  heating  efficiency  of  nanoparticle 
solutions. 
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Figure  7.  Temperature  traces  for 
three  heating/cooling  cycles  of  an 
aqueous  solution  of  5mM  6:1- 
aspect-ratio  nanorods  (red  curve) 
compared  to  pure  DI  water  (black 
curve).  The  sample  was  illuminated 
with  microwave  energy  for  three 
minutes  for  each  cycle  and  turned 
off  for  four  minutes  for  the  sample 
to  cool. 


Task  3.4c:  SWCNTs.  We  conducted  heating  efficiency  experiments  (at  3  GHz)  using  the  TM 
materials  described  in  Task  3.1,  mixed  with  various  concentrations  of  SWCNTs  that  have  been 
identified  to  be  non-toxic  in  mice.  We  chose  3  GHz  because  of  its  relevance  to  both  microwave- 
induced  thermoacoustic  imaging  and  hyperthermia  treatment.  We  prepared  the  samples  for  the 
heating  experiments  by  pouring  the  liquid  TM  mixture  into  a  1.1-mm-inner-diameter  glass 
capillary  tube.  The  mixture  was  allowed  to  gel  around  a  fiber-optic  temperature  probe  connected 
to  a  calibrated  fluoroptic  thermometer.  We  note  that  the  reference  sample  was  constructed  from 
a  pure  TM  mixture  with  the  same  amount  of  Pluronic  as  the  other  samples,  but  with  no 
SWCNTs.  The  capillary  tube  was  inserted  through  a  small  hole  drilled  into  the  center  of  the 
broad  wall  of  an  S-band  waveguide.  A  microwave  synthesizer  and  amplifier  generated  1  W  of 
continuous  microwave  power  at  3  GHz  that  was  delivered  to  the  sample  via  the  waveguide.  The 
fluoroptic  thermometer  recorded  the  time-varying  temperature  of  each  sample  as  it  was  heated. 
We  heated  each  sample  for  3  min,  then  turned  off  the  microwave  source  and  allowed  it  to  cool 
for  5  min.  The  complete  set  of  results  from  this  study  is  reported  in  Appendix  C  (Figure  3  and 
Table  II  therein).  We  demonstrated  measurement  repeatability  for  each  concentration.  The 
measured  data  shows  that  as  the  concentration  of  SWCNTs  increases,  the  maximum  temperature 
also  increases.  The  temperature  increase  scales  roughly  linearly  with  the  increase  in  conductivity 
reported  in  Task  3.1  (and  in  Table  I  of  Appendix  C).  SWCNT  concentrations  as  small  as  0.22% 
by  weight  led  to  an  average  steady-state  temperature  rise  that  was  6°C  higher  than  the  rise 
observed  in  the  TM  material  without  SWCNTs.  Similarly  promising  results  were  obtained  for  the 
dispersions  of  SWCNTs  in  water  and  PL-PEG  (Appendix  E). 
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Key  Research  Accomplishments 

•  Validated  two  complementary  techniques  for  microwave-frequency  dielectric  properties 
measurements,  namely  a  wideband  (0.5-20  GHz)  open-ended  coaxial  probe  technique  and  a 
narrowband  cavity  perturbation  technique. 

•  Conducted  wideband  and  narrowband  microwave-frequency  measurements  of  the  effective 
electromagnetic  properties  of  liquid  dispersions  of  candidate  micro/nanoparticle  contrast 
agents  including  air-filled  microbubbles,  metallic  nanoparticles,  and  carbon  nanotubes. 

•  Demonstrated  that  air-filled  microbubbles  with  clinically  relevant  concentrations 
significantly  decrease  the  microwave  dielectric  properties  of  the  liquid  dispersion  relative  to 
the  pure  liquid  without  microbubbles,  and  concluded  that  microbubbles  are  a  promising 
contrast  agent  for  contrast-enhanced  microwave  breast  cancer  detection. 

•  Demonstrated  that  compact  metallic  nanoparticles  do  not  significantly  alter  the  microwave 
dielectric  properties  of  the  liquid  dispersion  relative  to  the  pure  liquid  without  nanoparticles, 
and  therefore  are  unlikely  to  be  promising  as  contrast  agents  for  contrast-enhanced 
microwave  breast  cancer  detection. 

•  Demonstrated  that  single-wall  carbon  nanotubes  (SWCNTs)  significantly  increase  the 
microwave  dielectric  properties  of  the  liquid  dispersion  relative  to  pure  liquid  without 
SWCNTs.  Also  demonstrated  that  SWCNTs  at  biologically  relevant  concentrations 
significantly  increase  the  microwave  dielectric  properties  of  tumor-micking  phantom 
materials,  and  concluded  that  SWCNTs  are  a  promising  contrast  agent  for  contrast-enhanced 
microwave  breast  cancer  detection. 

•  Developed  2D  classical  electromagnetics  simulation  tools  for  calculating  the  effective 
dielectric  properties  of  micro/nanoparticle  dispersions  and  used  the  tools  to  investigate  a 
variety  of  micro/nanoparticle  dispersions  at  microwave  frequencies. 

•  Compared  2D  numerical  simulation  results  with  experimental  measurements  of  microbubble 
dispersions  and  obtained  excellent  agreement  in  the  computed  and  measured  effective 
dielectric  properties. 

•  Compared  2D  numerical  simulation  results  with  experimental  measurements  of  dispersions 
containing  compact  (isotropic  or  low-aspect-ratio  anisotropic)  metallic  nanoparticles  and 
found  discrepancies  in  the  computed  and  measured  values  that  may  be  attributed  to 
interfacial  effects  and/or  limitations  of  the  2D  simulation. 

•  Developed  a  2D  numerical  simulation  tool  that  couples  the  full-wave  FDTD  method  for 
solving  Maxwell’s  equations  with  an  ensemble  Monte  Carlo  approach  for  solving  the 
Boltzmann  transport  equation,  and  validated  the  tool  by  simulating  electromagnetic  wave 
interactions  with  low-doped  silicon  at  THz  frequencies. 

•  Characterized  the  dielectric  properties  of  a  variety  of  CNT-based  particles  in  liquid 
dispersions  and  identified  single -wall  CNTs  with  large  aspect  ratios  (e.g.  length-to-diameter 
ratios  on  the  order  of  -1000:1)  to  offer  the  greatest  field  enhancement  and  corresponding 
increase  in  effective  dielectric  properties. 
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•  Developed  a  fabrication  procedure  for  constructing  tumor-mimicking  phantoms  with  and 
without  SWCNT  mixtures  at  concentrations  known  to  be  non-toxic  in  mice. 

•  Developed  state-of-the-art  synthesis  protocols  of  spherical  nanoparticles  and  nanorods  and 
modified  their  surface  properties  to  stably  suspend  them  in  saline  solutions  representing 
realistic  physiological  conditions. 

•  Conducted  microwave  heating  experiments  on  low-concentration  dispersions  of  compact 
metallic  nanoparticles  and  characterized  their  microwave  heating  efficiency  relative  to  pure 
liquid  samples  without  nanoparticles.  Demonstrated  that  the  observed  increase  in  heating 
efficiency  was  due  to  the  surfactant  rather  than  the  nanoparticles  themselves. 

•  Conducted  microwave  heating  experiments  on  low-concentration  dispersions  of  SWCNTs 
and  characterized  their  microwave  heating  efficiency  relative  to  pure  liquid  samples  with  the 
surfactant  but  without  nanotubes.  Demonstrated  that  the  observed  increase  in  heating 
efficiency  was  indeed  due  to  the  CNTs. 

•  Demonstrated  that  microbubbles  diminish  the  thermoacoustic  response  of  simple  “tumor” 
targets  compared  to  the  pre-contrast-agent  response,  indicating  that  such  agents  may  enhance 
the  sensitivity  of  cancer  detection  through  the  use  of  differential  thermoacoustic  imaging 
techniques. 
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Conclusions 

During  the  course  of  this  project,  we  achieved  significant  milestones  in  our  contrast- 
enhanced  microwave  breast  cancer  detection  research  -  the  demonstration  of  significant  changes 
in  the  microwave-frequency  effective  dielectric  properties,  heating  efficiency,  and 
thermoacoustic  response  of  mixtures  containing  either  microbubbles  or  single-wall  carbon 
nanotubes.  To  the  best  of  our  knowledge,  these  accomplishments  represent  the  first 
experimental  laboratory  demonstrations  of  contrast  enhancement  in  a  context  relevant  to 
microwave  breast  cancer  detection.  Our  results  have  advanced  our  understanding  of  microwave 
interactions  with  microparticles  and  nanoparticles  in  dispersions  and  have  revealed  that  a  very 
large  aspect  ratio  is  required  in  order  for  metallic  particles  to  enhance  the  microwave  response. 

In  a  complementary  NIH-funded  research  project  involving  the  development  of  microwave 
breast  imaging  techniques,  we  have  investigated  the  impact  of  the  level  of  dielectric-properties 
changes  found  in  our  contrast  agent  research  (this  project)  on  the  sensitivity  of  detecting  small 
(~1  cm)  tumors  in  anatomically  realistic  numerical  breast  phantoms.  We  created  3D  microwave 
images  of  each  breast  phantom  by  applying  inverse  scattering  algorithms  to  the  microwave 
scattered  signals  recorded  for  each  phantom,  both  before  and  after  the  contrast  agent  injection. 
We  are  very  encouraged  by  the  fact  that  tumors  were  detected  via  differential  imaging  in 
phantoms  representating  a  wide  range  of  breast  tissue  densities,  ranging  from  mostly  fatty  to 
extremely  dense. 

Thus  our  work  to  date  strongly  suggests  that  microwave  imaging  enhanced  with 
microparticle  or  nanoparticle  contrast  agents  svstemically  delivered  to  the  tumor  has  the  potential 
to  overcome  many  of  the  limitations  of  conventional  breast  cancer  screening  modalities.  The 
physical  basis  for  our  microwave  imaging  approach  is  the  combined  intrinsic  and  extrinsic 
microwave-frequency  dielectric  contrast  between  normal  breast  tissue  and  malignant  tumors 
containing  metal  nanoparticles.  Our  proposed  approach  takes  advantage  of  biophysical  contrast 
mechanisms,  such  as  water  content,  angiogenesis,  blood-flow  rate,  and  temperature,  and 
dramatically  augments  them  by  the  selective  introduction  of  metallic  nanoparticles  into  the  tumor 
via  passive  or  targeted  delivery.  Microwave  imaging  techniques  that  exploit  the  resultant 
unprecedented  contrast  should  provide  the  sensitivity  and  resolution  needed  to  reliably  detect 
extremely  small  (less  than  1  cm)  malignant  tumors  even  in  challenging  cases  of  radiographically 
dense  breast  tissue  or  tumors  located  in  the  upper  outer  breast  quadrant  near  the  chest  wall  where 
an  estimated  50%  of  all  breast  tumors  occur.  Microwave  imaging  can  also  take  advantage  of 
enhanced  backscatter  due  to  vascularization  and  accumulation  of  contrast  agent  particles  in 
malignant  tumors,  as  well  as  morphology-dependent  characteristics  such  as  spectral  and 
polarization  signatures,  for  the  purpose  of  lesion  discrimination  and  characterization.  The 
predicted  safety,  comfort  (no  breast  compression),  and  ultra-low-cost  features  of  contrast- 
enhanced  microwave  imaging  means  that  this  breast  cancer  screening  technology  could  be  made 
widely  available  to  the  general  public,  including  medically  under-served  populations,  and  should 
improve  public  compliance  with  annual  screening  recommendations. 
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Abstract 

Microwave-induced  thermoacoustic  tomography  (MI-TAT)  is  an  imaging 
technique  that  exploits  dielectric  contrast  at  microwave  frequencies  while 
creating  images  with  ultrasound  resolution.  We  propose  the  use  of 
microbubbles  as  a  dielectric  contrast  agent  for  enhancing  the  sensitivity  of 
MI-TAT  for  breast  cancer  detection.  As  an  initial  investigation  of  this 
concept,  we  experimentally  studied  the  extent  to  which  the  microwave-induced 
thermoacoustic  response  of  a  dielectric  target  is  modified  by  the  presence  of  air- 
filled  glass  microbubbles.  We  created  mixtures  of  ethylene  glycol  with  varying 
weight  percentages  of  microbubbles  and  characterized  both  their  microwave 
properties  (0.5-6  GHz)  and  thermoacoustic  response  when  irradiated  with 
microwave  energy  at  3  GHz.  Our  data  show  that  the  microbubbles  considerably 
lowered  the  relative  permittivity,  electrical  conductivity  and  thermoacoustic 
response  of  the  ethylene  glycol  mixtures.  We  hypothesize  that  the  interstitial 
infusion  of  microbubbles  to  a  tumor  site  will  similarly  create  a  smaller 
thermoacoustic  response  compared  to  the  pre-contrast- agent  response,  thereby 
enhancing  sensitivity  through  the  use  of  differential  imaging  techniques. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Microwave-induced  thermoacoustic  tomography(MI-TAT)  is  a  hybrid  imaging  modality 
that  exploits  dielectric-properties  contrasts  while  creating  images  with  ultrasound-quality 
resolution.  In  thermoacoustic  imaging,  the  tissue  is  irradiated  with  sub-microsecond 
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electromagnetic  pulses.  This  electromagnetic  energy  is  selectively  absorbed  by  the  higher- 
conductivity  tissues  and  induces  acoustic  waves  by  the  means  of  thermoelastic  expansion. 
These  thermoacoustic  signals  are  detected  by  ultrasound  transducers  and  processed  for  image 
reconstruction.  One  promising  application  for  microwave  imaging  techniques  such  as  MI-TAT 
(as  well  as  others)  is  early-stage  breast  cancer  detection  (Kruger  et  al  1999,  2002,  Ku  et  al 
2005,  Xu  et  al  2005,  Meaney  et  al  2000,  Li  et  al  2005). 

A  recently  published  large-scale  spectroscopy  study  conducted  by  the  University  of 
Wisconsin  and  the  University  of  Calgary  (Lazebnik  et  al  2007)  has  shown  that  the  contrast  in 
the  microwave-frequency  dielectric  properties  between  malignant  tissues  and  healthy  adipose- 
dominated  tissues  in  the  breast  is  as  large  as  10:1,  but  the  contrast  between  malignant  and 
normal  glandular/fibroconnective  tissues  in  the  breast  is  not  more  than  about  10%.  Nearly  all 
breast  cancers  originate  in  the  ducts  or  lobules  of  the  breast,  which  are  glandular  tissues.  The 
use  of  contrast  agents,  such  as  those  that  accumulate  in  tumors  via  passive  mechanisms  (Maeda 
et  al  2000),  may  enhance  the  dielectric  contrast  between  malignant  and  normal  glandular  tissue, 
thereby  increasing  the  sensitivity  of  microwave  imaging  for  breast  cancer  detection. 

Here,  we  propose  air- filled  microbubbles  as  a  contrast  agent  in  MI-TAT.  Because  of 
their  echo-enhancing  properties,  microbubbles  have  been  studied  for  the  past  several  decades 
as  a  contrast  agent  in  ultrasound  imaging  (Calliada  et  al  1998,  Cosgrove  2006,  Wells 
2006).  Microbubbles  may  also  be  advantageous  in  MI-TAT  modalities  because  of  their  low 
microwave  absorption.  We  hypothesize  that  the  thermoacoustic  response  of  a  tumor  infused 
with  microbubbles  will  be  not  only  smaller  than  that  of  surrounding  normal  glandular  tissue 
structures  but  also  smaller  than  the  response  of  the  tumor  in  its  pre-contrast-agent  state.  This 
extrinsic  dielectric  contrast  mechanism  could  be  exploited  via  differential  imaging  whereby 
the  breast  is  imaged  before  and  after  the  contrast  agent  is  injected. 

As  an  initial  step  toward  evaluating  the  feasibility  of  microbubbles  as  a  contrast  agent  in 
MI-TAT,  we  tested  the  extent  to  which  the  microwave-induced  response  of  a  simple  dielectric 
target  is  modified  by  the  presence  of  microbubbles.  The  remainder  of  this  paper  is  organized  as 
follows.  Section  2  summarizes  the  theory  of  thermoacoustic  phenomena.  Section  3  describes 
the  methods  used  to  create  the  thermoacoustic  targets  and  experimentally  characterize  both 
their  dielectric  properties  and  thermoacoustic  response.  Section  4  discusses  the  results  of  the 
study,  and  the  conclusions  are  summarized  in  section  5. 


2.  Theory 

Thermoacoustic  theory  has  been  detailed  in  a  number  of  review  papers  (see,  for  example,  Tam 
1986  and  Xu  and  Wang  2006).  Here,  we  briefly  summarize  the  fundamental  concepts.  The 
thermoacoustic  differential  equation  given  by  equation  (1)  is  derived  from  the  equations  of 
motion  and  thermal  expansion  assuming  a  homogeneous  acoustic  and  thermal  medium 

V2p(r,t)  -  U?(r,0  =  -pep^T(r,t).  (1) 

Here  p(r,t )  is  the  acoustic  pressure  at  position  r  and  time  t,  /3e  is  the  volume  expansion 
coefficient,  T(r,t )  is  the  temperature  rise  throughout  the  medium  due  to  microwave 
absorption,  p  is  the  mass  density  and  c  is  the  speed  of  sound. 

For  the  case  in  which  microwave  power  dissipation  in  a  region  is  uniform  and  so  rapid 
that  thermal  diffusion  can  be  neglected,  the  temperature  rise  in  that  region  can  be  described  as 

SAR-  r 
Cp 


AT  = 


(2) 
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Here  r  is  the  microwave  pulse  length,  Cp  is  the  specific  heat  and  S  AR  is  the  specific  absorption 
rate,  defined  as  a  \  E  |2/2p,  where  E  is  the  electric  field,  a  is  the  effective  electrical  conductivity 
and  p  is  the  mass  density. 

These  equations  suggest  that  a  lesion  with  lowered  conductivity  should  have  a  smaller 
local  temperature  rise  and  a  smaller  thermoacoustic  response. 


3.  Materials  and  methods 

3.1.  Thermoacoustic  targets 

Our  thermoacoustic  targets  were  composed  of  mixtures  of  ethylene  glycol  (EG)  with  0%, 
5%,  10%,  20%  and  30%  by  weight  concentrations  of  microbubbles.  These  concentrations  are 
within  the  typical  109  microbubbles/mL  dosage  administered  in  human  ultrasound  procedures 
(Stride  and  Saffari  2005).  The  microbubbles  used  were  air-filled  glass  spheres  (iM30K,  3M) 
with  an  average  diameter  of  18  jam  and  a  density  of  0.6  g  mL-1 .  Ethylene  glycol  was  chosen 
because  it  is  a  relatively  inexpensive  and  readily  available  liquid  with  dielectric  properties  that 
are  approximately  within  the  range  of  biological  tissues  at  the  microwave  frequency  of  interest 
(3  GHz).  In  this  study,  we  were  concerned  only  with  the  role  of  the  dielectric  properties  on 
the  thermoacoustic  response  and  did  not  consider  or  characterize  the  effect  of  the  thermal, 
mechanical  or  acoustic  properties  of  the  mixtures  or  their  constituents.  Thus  those  other 
properties  did  not  influence  our  choice  of  the  background  liquid. 

The  mixtures  were  prepared  by  vigorously  mixing  EG  with  the  microbubbles  for  2  min 
before  taking  any  measurements.  The  thermoacoustic  targets  were  constructed  by  pouring  the 
mixtures  into  plastic  tubes  with  an  inner  diameter  of  6  mm.  The  use  of  a  plastic  tube  for  target 
construction  has  been  used  in  experimental  work  reported  by  other  groups  (Ku  et  al  2005,  Jin 
et  al  2007,  Paltauf  et  al  2007).  Both  ends  of  the  tube  were  sealed  with  vinyl  plastic  putty.  The 
tube  was  placed  in  a  stand  and  immersed  in  an  oil  tank. 


3.2.  Dielectric  characterization  of  the  targets 

The  dielectric  properties  of  the  mixtures  were  characterized  using  a  well  established  open- 
ended  coaxial  probe  technique  described  in  Popovic  et  al  (2005).  Briefly,  the  tip  of  hermetically 
sealed,  stainless-steel/borosilicate  glass  coaxial  probe  was  immersed  in  the  liquid  under  test. 
The  complex  reflection  coefficient  was  recorded  using  a  performance  network  analyzer  and 
converted  into  the  complex  dielectric  properties  using  the  procedure  described  in  Popovic 
et  al  (2005).  As  reported  in  Popovic  et  al  (2005),  the  measurement  uncertainty  of  this 
technique  is  no  more  than  approximately  10%. 

We  conducted  dielectric  spectroscopy  measurements  on  freshly  mixed  EG/microbubble 
mixtures  contained  within  glass  test  tubes  of  the  same  height  as  the  plastic  tubes  used  in  the 
thermoacoustic  measurements.  Five  dielectric  measurements  were  made  at  three  different 
depths  in  the  test  tubes  spanning  the  frequency  range  0.5-6  GHz  and  the  results  at  each 
frequency  were  averaged.  The  range  of  depths  spanned  the  region  of  the  target  that  was 
illuminated  with  the  microwaves  in  the  thermoacoustic  experiments.  The  measurements  were 
made  at  room  temperature  (approximately  22  °C). 

Since  the  air-filled  microbubbles  are  less  dense  than  the  EG,  they  separate  from  the  mixture 
over  time  and  rise  to  the  top.  To  ensure  that  the  state  of  the  mixture  was  the  same  for  both  our 
dielectric  measurements  and  thermoacoustic  experiments,  we  made  dielectric  measurements 
2.5  min  after  mixing  the  solution — the  amount  of  time  that  elapsed  in  the  subsequent 
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RF  Enclosure 


Figure  1.  Thermoacoustic  experimental  setup. 


thermoacoustic  experiments  between  mixing  the  solution  and  recording  the  thermoacoustic 
response. 

3.3.  Thermoacoustic  characterization  of  the  targets 

A  schematic  of  the  thermoacoustic  experimental  setup  is  shown  in  figure  1 .  Microwave  pulses 
with  a  peak  power  of  approximately  30  kW  and  a  frequency  of  3  GHz  were  delivered  to  the 
thermoacoustic  target  via  an  S-band  rectangular  waveguide  (cross-sectional  dimensions  of 
72  mm  x  34  mm).  A  pulse  generator  (BNC  565,  Berkeley  Nucleonics)  was  used  to  set  the 
pulse  width  to  0.9  /xs  with  a  repetition  frequency  of  900  Hz.  These  settings  correspond  to  an 
average  power  of  24  W.  The  ultrasound  transducer  was  an  unfocused  piezoelectric  transducer 
with  a  2.25  MHz  center  frequency  (60%  6  dB  bandwidth)  and  12.7  mm  diameter  (V306, 
Olympus).  The  transducer  was  fastened  to  a  xyz-translation  stage  and  was  aligned  to  the 
target  using  the  pulse  echo  mode  of  the  ultrasound  pulser/receiver  (5800PR,  Olympus).  The 
ultrasound  pulser/receiver  was  shielded  from  the  microwave  generator  by  an  RF  enclosure  to 
eliminate  microwave  pickup. 

Safflower  oil  was  used  as  the  acoustic  coupling  medium.  The  use  of  oil  is  advantageous 
because  it  has  both  low  relative  permittivity  and  conductivity  and  therefore  provides  a  good 
dielectric  contrast  with  the  thermoacoustic  target  and  good  penetration  depth  into  the  tank. 
The  stand  and  walls  of  the  tank  were  constructed  out  of  acrylic,  which  has  similar  dielectric 
properties  as  oil. 

The  thermoacoustic  signal  detected  by  the  transducer  during  an  experimental  trial  with  a 
specific  target  was  amplified  by  60  dB  using  the  pulser/receiver  and  averaged  200  times  by 
the  oscilloscope  (Infiniium  54854A,  Agilent).  Both  the  microwave  generator  and  oscilloscope 
were  triggered  by  the  pulse  generator.  The  transducer  and  stand  were  kept  stationary  between 
different  thermoacoustic  measurements. 

3.4.  Speed  of  sound  characterization  of  the  targets 

We  made  speed  of  sound  measurements  using  the  well  accepted  substitution  technique 
described  in  Kremkau  et  al  (1981)  at  2.25  MHz.  Mixtures  of  EG  with  specific  microbubble 
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Figure  2.  (a)  Relative  permittivity  and  (b)  effective  conductivity  of  EG  solutions  with  varying 
concentrations  of  microbubbles  along  with  the  properties  of  non-infiltrated  fat  and  muscle  described 
by  the  four-pole  Cole-Cole  model  (Gabriel  et  al  1996). 


concentrations  were  poured  into  a  cylindrical  acrylic  tube  with  a  7.6  cm  diameter  and  a  25  fim 
polyvinylidene  chloride  (Saran  Wrap)  window.  The  sample  was  then  placed  in  a  water  tank 
and  illuminated  with  an  ultrasound  beam.  The  propagation  time  was  recorded  with  and  without 
the  sample.  The  speed  of  sound  was  then  calculated  from  the  relation  (Kremkau  et  al  1981): 


L.  -  .  \  U ) 

d  —  cwAt 

Here  d  is  the  length  of  the  cylindrical  sample,  cw  is  the  speed  of  sound  in  water  at  22  °C 
and  At  is  the  change  in  propagation  time  with  and  without  the  sample.  The  speed  of  sound 
measurements  were  made  2.5  min  after  mixing  the  solutions,  for  the  reasons  described  in 
section  3.2. 

4.  Results  and  discussion 

4.1.  Dielectric  properties 

Figures  2(a)  and  (b)  show  the  relative  permittivity  and  effective  conductivity  of  EG  and  the 
four  different  microbubble  solutions  over  a  frequency  range  of  0.5-6  GHz  at  the  time  that 
corresponded  to  when  the  thermoacoustic  response  was  recorded.  The  curves  depict  the 
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Table  1.  Summary  of  data  presented  in  figure  2. 


Initial 

microbubble 

concentration 
Solution  (%  by  weight) 

Estimated 

actual 

concentration 
(%  by  weight) 

Average 
cr  (3  GHz) 

Average 
or  (3  GHz) 
(S/m) 

1 

0 

0 

14.03 

2.32 

2 

5 

20 

9.66 

1.33 

3 

10 

30 

8.34 

1.11 

4 

20 

35 

7.38 

0.95 

5 

30 

40 

6.86 

0.84 

6 

100 

100 

NAa 

NAa 

a  NA:  not  available. 


average  of  the  different  measurements  and  the  vertical  bars  represent  the  deviation  of  the 
individual  measurements  from  the  average  value.  Table  1  summarizes  the  relative  permittivity 
and  conductivity  of  these  solutions  at  3  GHz.  Figure  2  shows  that  the  relative  permittivity  and 
conductivity  of  the  solutions  decrease  with  increasing  microbubble  concentration,  as  expected. 
In  addition  to  the  microwave  properties  of  the  solutions  investigated,  the  microwave  properties 
of  non-infiltrated  fat  and  muscle  are  plotted  using  the  four-pole  Cole-Cole  models  reported 
by  Gabriel  et  al  (1996)  for  comparison. 

Since  the  microbubbles  tend  to  float  to  the  top  of  the  mixture  over  time,  the  actual 
concentration  of  microbubbles  seen  by  the  microwaves  is  higher  than  that  of  the  original 
mixture.  In  order  to  estimate  the  actual  concentration  of  microbubbles  at  the  time  the 
thermoacoustic  response  was  recorded,  additional  higher  concentration  microbubble  mixtures 
were  made.  The  dielectric  properties  of  these  additional  mixtures  were  measured  immediately 
after  the  solution  was  created.  These  instantaneous  dielectric  measurements  of  the  higher 
concentration  mixtures  (not  shown)  were  then  compared  to  those  of  figure  2  to  determine  the 
actual  weight  concentration  of  the  mixtures  at  the  time  of  the  thermoacoustic  measurements. 
The  estimated  concentrations  are  summarized  in  table  1 .  These  measurements  establish  that 
the  dielectric  properties  of  a  solution  that  was  created  by  using  5%  by  weight  concentration 
of  microbubbles,  at  two  and  a  half  minutes,  are  equivalent  to  a  20%  by  weight  concentration 
of  microbubbles.  Similarly,  mixtures  created  by  using  10%,  20%  and  30%  by  weight 
concentration  of  microbubbles,  at  two  and  a  half  minutes,  have  dielectric  properties  of  solutions 
with  30%,  35%  and  40%  by  weight  microbubbles,  respectively. 


4.2.  Thermoacoustic  response 

Figure  3  shows  the  temporal  thermoacoustic  waveform  as  detected  by  the  ultrasound  transducer 
for  Solution  1  (pure  EG  solution),  for  five  different  trials,  where  each  trial  involved  a  fresh 
mixture  of  the  same  concentration.  The  microwave  pulse  illuminates  the  target  at  t  =  0  /xs. 
Since  the  electromagnetic  energy  is  delivered  to  the  target  at  the  speed  of  light,  the  acoustic 
waves  are  induced  at  t  ~  0  /xs.  The  approximately  34  /xs  delay  corresponds  to  the  propagation 
time  of  the  induced  acoustic  waves  to  reach  the  transducer  in  oil.  Each  thermoacoustic  temporal 
waveform  shown  in  figure  3  is  composed  of  mainly  three  lobes.  This  is  the  expected  pressure 
profile  emitted  from  a  cylindrical  thermoacoustic  target  inside  a  background  medium  with  a 
different  density  and  speed  of  sound  (Diebold  et  al  1990).  The  first  positive  and  first  negative 
lobes  are  characteristic  of  two-dimensional  waves  emitted  from  a  cylindrical  acoustic  source. 


Microbubbles  for  enhancing  thermocoustic  contrast 


647 


Figure  3.  Thermoacoustic  temporal  waveform  of  Solution  1  for  five  independent  trials.  The 
waveform  shown  for  each  trial  represents  an  average  of  200  measurements.  The  small  variations 
observed  across  waveforms  reflect  slight  differences  between  trials  (e.g.,  small  differences  in  the 
plastic  tubes  and  their  positions). 


The  third  positive  lobe  and  ensuing  acoustic  signature  trail  is  due  to  the  difference  in  the 
density  and  speed  of  sound  of  the  target,  EG,  and  the  background  medium,  safflower  oil. 

We  used  the  multiple  trials  reported  in  figure  3  to  characterize  the  repeatability  of  our 
measurements  under  experimental  conditions  that  are  intended  to  be  invariant  (e.g.,  fixed 
concentration  of  microbubbles).  This  assessment  is  critical  for  correctly  interpreting  the 
differences  observed  in  the  thermoacoustic  waveform  under  intentionally  varying  experimental 
conditions  (differing  concentrations,  for  example).  Figure  3  shows  that  the  variations  in  the 
thermoacoustic  response  of  similar  dielectric  targets  are  no  greater  than  approximately  6% 
in  the  first  peak  and  10%  in  the  temporal  width  when  compared  to  the  average  values.  The 
temporal  fluctuations  are  no  greater  than  approximately  1  /xs.  We  attribute  these  variations  in 
the  thermoacoustic  waveform  to  small  changes  in  the  physical  properties  of  the  plastic  tube 
and  the  placement  of  the  target  relative  to  the  ultrasound  transducer.  These  variations  are 
inevitably  introduced  between  different  trials. 

Figure  4  shows  thermoacoustic  temporal  waveforms  for  one  representative  trial  of  each  of 
the  solutions.  The  average  and  the  range  of  minimum  and  maximum  values  of  the  first  peak  as 
well  as  the  temporal  width  of  the  thermoacoustic  response  for  the  five  different  trials  of  each 
concentration  are  summarized  in  table  2.  Solution  1  (pure  EG)  and  Solution  6  (the  plastic 
tube  filled  with  only  microbubbles)  serve  as  an  upper  and  a  lower  bound  of  the  thermoacoustic 
response.  The  thermoacoustic  response  of  Solution  6  is  very  small,  and  most  likely  due  to  the 
microwave  absorption  of  the  glass  shells  of  the  microbubbles  and  the  plastic  tube.  The  change 
observed  across  solutions  is  clearly  larger  than  the  variation  for  repeat  measurements  with  the 
same  solution,  as  seen  in  figure  3. 
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Figure  4.  Thermoacoustic  temporal  waveform  for  all  solutions. 


Table  2.  Summary  of  data  presented  in  figure  4  and  speed  of  sound  measurement  results. 


Average 
amplitude 
of  first  peak 
in  TAa 

waveform 
Solution  (a.u.) 

Min/max 
amplitude 
of  first  peak 
in  TAa 

waveform 

(a.u.) 

Average 
temporal 
widthb 
of  TAa 

waveform 

(/xs) 

Speed 
of  sound 
(m/s) 

1 

0.214 

0.203/0.227 

5.1 

1660.7 

2 

0.182 

0.171/0.202 

3.56 

1729.8 

3 

0.176 

0.164/0.188 

3.50 

1805.2 

4 

0.149 

0.138/0.155 

3.28 

1922.6 

5 

0.116 

0.103/0.127 

3.18 

2049.4 

6 

^0 

^0/^0 

2.6 

NAC 

a  TA:  Thermoacoustic. 

b  Measured  by  calculating  the  temporal  distance  between  the  two  positive  peaks  in  the 
thermoacoustic  waveform. 
c  NA:  not  available. 


Figure  4  and  the  data  in  table  2  show  that  the  thermoacoustic  response  diminishes  with 
increasing  concentrations  of  microbubbles.  Not  only  does  the  addition  of  air-microbubbles 
modify  the  dielectric  properties  (see  figure  2),  but  it  is  also  well  known  that  gas-filled 
microbubbles  increase  acoustic  attenuation  due  to  acoustic  scattering  effects  (Chen  et  al 
2002).  Therefore,  it  is  important  to  distinguish  whether  the  reduction  in  the  initial  peak  of 
the  measured  thermoacoustic  waveform  after  the  addition  of  microbubbles  is  due  to  acoustic 
attenuation  or  dielectric  properties  modification.  The  analytic  solution  for  the  thermoacoustic 
impulse  response  of  a  target  without  acoustic  attenuation  has  been  discussed  in  great  detail 
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in  Wang  and  Wu  (2007).  This  derivation  shows  that  the  leading  edge  of  the  thermoacoustic 
waveform  is  due  to  a  pressure  wave  emitted  from  the  target  surface  closest  to  the  detector. 
This  physical  fact  remains  true  even  if  the  non- attenuating  target  medium  is  replaced  with  an 
acoustically  lossy  target.  The  amplitude  of  the  acoustic  disturbance  from  the  closest  target 
surface  will  not  be  affected  by  any  acoustic  attenuation  within  the  target  volume.  Therefore, 
one  can  conclude  that  the  observed  reduction  in  the  first  peak  of  the  thermoacoustic  response 
with  increasing  concentrations  of  microbubbles  is  due  to  the  reduced  effective  conductivity 
of  the  mixtures  as  shown  by  the  dielectric  measurements.  The  lower  conductivity  leads  to 
smaller  microwave  dissipation  and  hence  to  a  lower  thermoacoustic  response. 

4.3.  Speed  of  sound  measurements 

Another  observation  from  figure  4  and  the  data  in  table  2  is  the  reduction  in  temporal  width 
of  the  thermoacoustic  waveforms  with  higher  concentrations  of  microbubbles.  The  temporal 
width  of  the  induced  thermoacoustic  wave  is  directly  related  to  the  acoustical  width  of  the 
target  (Diebold  et  al  1991).  As  discussed  in  the  review  article  by  Goldberg  et  al  (1993), 
several  research  groups  have  shown  that  the  speed  of  sound  of  a  medium  decreases  with 
increasing  concentrations  of  gas-filled  organic-shell  microbubbles.  Therefore,  the  expected 
trend  of  introducing  gas  microbubbles  is  that  the  thermoacoustic  pulse  width  should  widen 
since  the  speed  of  sound  in  air  is  much  slower  than  the  speed  of  sound  in  EG.  However,  in  our 
experiments,  we  observed  a  decrease  in  the  temporal  width  of  the  thermoacoustic  waveform 
with  increasing  concentrations  of  microbubbles,  suggesting  an  increase  in  the  speed  of  sound 
in  our  microbubble  mixtures. 

This  inference  was  confirmed  with  the  speed  of  sound  measurements  summarized  in 
table  2.  These  data  show  that  as  the  concentration  of  microbubbles  increases  the  speed 
of  sound  in  the  mixture  increases  as  well.  This  is  consistent  with  the  trends  observed  in 
figure  4  and  table  2  in  which  the  thermoacoustic  temporal  width  decreased  with  increasing 
concentration  of  microbubbles. 

5.  Conclusions 

In  these  experiments  we  characterized  the  dielectric  and  acoustic  properties  as  well  as 
the  thermoacoustic  response  of  mixtures  with  varying  concentration  of  air-filled  glass 
microbubbles.  Microbubbles  significantly  lower  the  microwave  absorption  of  the  target  which 
should  reduce  the  thermoacoustic  response.  These  microbubbles  also  increase  the  acoustic 
velocity  which  should  reduce  the  temporal  width  of  the  thermoacoustic  response.  In  our 
thermoacoustic  experiments  we  observed  reductions  of  both  the  magnitude  and  temporal 
width  of  the  thermoacoustic  response  with  higher  levels  of  microbubble  concentrations.  With 
organic-shell  microbubbles  one  can  expect  the  same  reduction  in  amplitude,  but  a  broadening 
of  the  temporal  width.  The  changes  in  amplitude  of  the  thermoacoustic  response  are  promising 
enough  to  warrant  further  investigation  into  the  use  of  microbubbles  as  a  potential  contrast 
agent  in  MI-TAT. 
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Abstract  The  interactions  between  carriers  and  fields  in 
semiconductors  at  low  frequencies  (<100  GHz)  can  be  ad¬ 
equately  described  by  numerical  solution  of  the  Boltzmann 
transport  equation  coupled  with  Poisson’s  equation.  As  the 
frequency  approaches  the  THz  regime,  the  quasi-static  ap¬ 
proximation  fails  and  full-wave  dynamics  must  be  consid¬ 
ered.  Here,  we  review  recent  advances  in  global  modeling 
techniques — numerical  techniques  that  couple  carrier  dy¬ 
namics  with  full  wave  dynamics.  We  focus  on  the  coupling 
between  the  stochastic  ensemble  Monte  Carlo  (EMC)  sim¬ 
ulation  of  carrier  transport  and  the  finite-difference  time- 
domain  (FDTD)  solution  to  Maxwell’s  curl  equations.  We 
discuss  the  stability  and  accuracy  requirements  for  different 
types  of  high-frequency  excitation  (wave  illumination  vs. 
ac  bias),  and  present  simulation  results  for  the  THz-regime 
conductivity  of  doped  bulk  silicon,  ultrafast  carrier  dynam¬ 
ics  and  radiation  patterns  in  GaAs  filaments,  and  the  ac  re¬ 
sponse  of  GaAs  MESFETs. 

Keywords  Global  modeling  •  THz  conductivity  • 
Microwave  devices  •  EMC-FDTD  •  Monte  Carlo  simulation 


1  Introduction 

Typical  doped  semiconductors  have  plasma  frequencies 
and  characteristic  carrier  scattering  rates  in  the  high  mi- 
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crowave  and  terahertz  (THz)  frequency  regimes.  Outside  of 
the  regime  where  the  stimulating  frequency  is  comparable 
to  the  carrier  relaxation  rate,  the  effects  of  carrier  trans¬ 
port  and  electromagnetic  (EM)  wave  propagation  are  suf¬ 
ficiently  decoupled  to  permit  independent  treatment  in  nu¬ 
merical  models.  Accordingly,  most  of  the  recent  advances 
and  applications  of  charge  transport  solvers  use  simple  elec¬ 
trostatic  solvers  to  incorporate  electric  field  effects  [1-4]. 
At  the  same  time,  EM  analyses  use  simpler  bulk  material 
models  alongside  full- wave  electrodynamic  solvers  [5]. 

As  device  speeds  continue  to  increase  in  step  with  de¬ 
creasing  critical  dimensions,  electrodynamic  effects  directly 
influence  high-frequency  device  performance,  and  analyses 
that  rely  on  quasi-static  electric  fields  lose  accuracy  [6,  7]. 
Similarly,  as  the  frequency  of  stimulating  EM  fields  in¬ 
creases  from  the  lower  microwave  range,  the  assumptions 
inherent  to  low-frequency  material  models  (namely,  that 
wr  <  1,  where  co  is  the  stimulating  frequency  and  t-1  is  the 
scattering  rate)  lose  validity.  This  prominent  interplay  be¬ 
tween  carrier  dynamics  and  electromagnetic  wave  dynamics 
in  the  sub-THz  and  THz  regimes  has  been  the  driving  force 
behind  advancements  in  modeling  techniques  that  combine 
a  charge  transport  kernel  with  a  computational  electromag¬ 
netic  (CEM)  solver,  hereafter  referred  to  as  global  modeling 
techniques. 

Interest  in  plasma  fusion  devices  prompted  early  work 
into  the  development  of  global  models.  Pioneering  work  by 
Langdon  and  Dawson  demonstrated  for  the  first  time  the  vi¬ 
ability  of  an  EM-particle  solver  [8].  Buneman’s  1968  report 
[9]  detailed  the  advantages  of  time-integrating  Maxwell’s 
equations  using  the  Yee  cell  [10].  Several  applications  of  this 
technique  followed  [11-13].  Boris  gave  a  detailed  account 
of  the  global  model  as  it  stood  in  1970  [14].  Advances  in  the 
1970s  addressed  self-force  [15]  and  numerical  Cherenkov 


Published  online:  12  August  2009 


4?)  Springer 


J  Comput  Electron 


instabilities  [16,  17].  The  global  models  resulting  from  these 
and  other  advances  are  described  in  detail  in  Refs.  [18-20]. 

The  technique  was  limited  by  the  computational  capabil¬ 
ities  of  the  era.  EM  calculations  exhibited  a  continual  linear 
increase  in  noise  [21,  22].  Additionally,  system  timescales 
varied  too  greatly  to  include  major  relevant  interactions 
without  forcing  a  reduced  speed  of  light  and  simultaneously 
increased  electron  mass.  A  1980  review  by  Buneman  et  al. 
suggested  that  a  3D  EM-particle  tracking  solver  could  pro¬ 
vide  qualitative  insight,  but  quantitative  results  were  only 
possible  for  reduced  systems  [23]. 

For  a  number  of  years  after  the  early  1980s,  parti¬ 
cle  tracking  models  and  CEM  solvers  advanced  more  or 
less  independently.  Many  of  the  early  EM-particle  re¬ 
ports  used  electromagnetic  models  that  were  clearly  early 
finite-difference  time-domain  (FDTD)  implementations. 
Researchers  continued  to  advance  FDTD;  a  detailed  history 
is  available  in  Refs.  [5,  24].  The  Monte  Carlo  technique  had 
already  been  of  great  interest  for  many  years,  and  the  disci¬ 
pline  continued  to  mature  in  this  time;  for  a  history  of  the 
ensemble  Monte  Carlo  (EMC)  technique  see  Ref.  [25]. 

The  need  for  greater  efficiency  and  decreased  memory  re¬ 
quirements  drove  exploration  into  alternative  global  model¬ 
ing  approaches.  Hydrodynamics  transport  models  and  drift- 
diffusion  models  have  been  used  successfully  in  place  of 
EMC  in  many  cases  [7,  26-35]  offering  decreased  com¬ 
putational  overhead  [7].  Numerical  simulations  combining 
full-wave  electromagnetic  solutions  via  modern  FDTD  with 
particle-based  transport  models  via  EMC  were  first  reported 
by  El-Ghazaly,  Joshi,  and  Grondin  in  1990  [36].  This  work 
provided  a  more  accurate  model  of  sub-picosecond  carrier 
transport  in  photoconductive  switches  than  had  previously 
been  available.  A  comprehensive  review  of  the  global  mod¬ 
eling  efforts  of  the  1980s  and  1990s  was  conducted  by 
Grondin,  El-Ghazaly,  and  Goodnick  [7] . 

Modern  improvements  in  computer  processing  power 
have  eased  the  restrictions  on  computational  burden  in  cur¬ 
rent  device  and  materials  simulations,  allowing  researchers 
to  use  more  computationally  intensive  techniques  such  as 
EMC  and  FDTD.  Application  of  a  combined  3D  full-band 
cellular  Monte  Carlo  (CMC)  device  simulator  with  a  3D 
full- wave  FDTD  solver  was  first  demonstrated  in  2003,  for 
modeling  the  electromagnetic  environment  surrounding  a 
simple  pin  diode  [37].  This  simulator  was  used  to  analyze 
high-frequency  transistor  behavior  using  direct  capture  of 
incident  and  reflected  voltage  waves  from  a  full-wave  analy¬ 
sis  of  a  3D  MESFET  device  [38].  Further  developments  in¬ 
clude  accurate  models  of  electron  transport  [39,  40]  and  the 
full-wave  effects  [41]  in  ultrasubmicrometer-gate,  pseudo- 
morphic,  high-electron  mobility  transistors  (pHEMTs).  Re¬ 
cently,  a  2D  combined  EMC-FDTD  solver  was  developed 
and  employed  for  high-frequency  characterization  of  doped 
silicon  [42,  43]. 


In  this  paper,  we  describe  recent  developments  in  the  field 
of  global  modeling  of  semiconductor  materials  and  devices, 
with  a  focus  on  coupled  EMC-FDTD  simulations.  In  Sect.  2, 
we  highlight  the  salient  features  of  EMC  and  FDTD  tech¬ 
niques,  with  the  assumption  that  the  reader  is  not  familiar 
with  both  techniques.  This  background  overview  is  followed 
by  details  on  the  implementation  of  the  combined  EMC- 
FDTD,  with  focus  on  the  importance  of  current  conserva¬ 
tion  and  the  definition  of  charge  assignment  and  current 
assignment  schemes.  The  application  of  this  global  solver 
to  doped  bulk  silicon  characterization  at  THz  frequencies 
is  presented  in  Sect.  3.  Excellent  agreement  between  the 
numerically-extracted  complex  conductivity  of  doped  sili¬ 
con  and  available  experimental  data  demonstrates  the  va¬ 
lidity  of  the  global  technique  for  high-frequency  semicon¬ 
ductor  analysis.  Section  4  describes  the  simulation  of  mi¬ 
crowave  transistors  using  the  global  CMC-FDTD  model. 
That  work  uses  implementation  modifications  necessary  for 
ac  device  analysis,  and  uses  a  faithful  description  of  ultrafast 
carrier  dynamics  and  radiation  patterns  to  characterize  thin 
intrinsic  GaAs  filaments  and  GaAs  MESFETs.  Finally,  in 
Sect.  5  we  summarize  these  results  and  provide  an  outlook 
on  the  future  of  global  modeling. 


2  Coupled  EMC-FDTD  simulation 

2.1  Ensemble  Monte  Carlo  (EMC)  simulation 


Ensemble  Monte  Carlo  (EMC)  is  a  powerful  stochastic 
method  used  for  numerical  simulation  of  carrier  trans¬ 
port  in  semiconductors  in  the  scattering-limited  (diffusive) 
regime  [25].  It  has  been  used  for  almost  four  decades  to  ac¬ 
curately  simulate  carrier  transport  properties  of  bulk  semi¬ 
conductors  and  semiconductor-based  devices,  and  it  pro¬ 
vides  a  benchmark  for  drift-diffusion  and  hydrodynamics 
equations  approaches  [2].  EMC  technique  details  are  dis¬ 
cussed  below;  for  more  information  see  Refs.  [1-4,  18,  25, 
44,  45]. 

On  timescales  longer  than  typical  relaxation  times,  EMC 
yields  the  numerical  solution  to  the  Boltzmann  transport 
equation  (BTE) 


—  +  v  ■  V?/  +  F  ■  V^/  =  — 
dt  rJ  pJ  dt 


(1) 


where  /(r,  p ,  t)  is  the  semiclassical  distribution  function,  v 
is  the  carrier  velocity,  and  F  is  the  total  force  acting  on  the 
carrier.  The  distribution  function  describes  the  instantaneous 
probability  that  a  carrier  exists  with  position  r  and  momen¬ 
tum  p.  f(r,p,t )  evolves  in  time  according  to  (1)  as  a  re¬ 
sult  of  diffusive  processes  (included  in  the  second  term)  and 
carrier  drift  due  to  external  forces  (described  by  the  third 
term).  The  collision  integral  (described  by  the  term  on  the 
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Fig.  1  Flowchart  of  the 
ensemble  Monte  Carlo  transport 
kernel 
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right  hand  side  of  the  equation)  describes  the  impact  of  the 
material- specific  scattering  mechanisms  on  the  carrier  en¬ 
semble. 

The  EMC  simulates  carrier  dynamics  in  semiconductors 
by  tracking  the  evolution  of  a  large  ensemble  of  particles 
[typically  O(105)]  through  time.  Each  carrier  undergoes  a 
series  of  scattering  events  and  free  flights.  Free-flight  mo¬ 
mentum  is  updated  according  to  the  Lorentz  force, 

„  dv  -> 

m  — =q(E  +  vxB),  (2) 

d  t 

where  m*  is  the  carrier  effective  mass,  q  the  carrier  charge, 

v  the  carrier  velocity,  and  E  and  B  are  the  electric  and 
magnetic  fields,  respectively.  A  random  number  generator  is 
used  to  calculate  the  duration  of  each  free  flight,  choose  the 
mechanism  for  the  next  scattering  event,  and  update  the  par¬ 
ticle’s  momentum  and  energy  state  as  needed,  according  to 
the  appropriate  statistical  probabilities.  Macroscopic  quan¬ 
tities  of  interest  (such  as  charge  density  and  drift  velocity) 
may  be  readily  extracted  via  ensemble  averages. 

The  majority  of  EMC  implementations  employ  the  quasi¬ 
static  assumption,  where  the  electromagnetic  period  of  os¬ 
cillation  is  sufficiently  long  compared  to  carrier  scattering 


times  that  the  field  may  be  considered  constant  within  any 
time  step.  A  simplified  flowchart  of  the  general  electrosta¬ 
tic  EMC  approach,  self-consistently  coupled  to  a  Poisson’s 
equation  solver,  is  shown  in  Fig.  1.  In  this  implementation, 
the  Poisson  solver  calculates  the  grid-based  scalar  potential 
that  results  from  the  instantaneous  charge  density  and  im¬ 
posed  biasing  conditions.  Poisson  solver  accuracy  requires 
a  grid  spacing  small  enough  to  resolve  the  smallest  feature 
of  interest  [46],  typically  chosen  to  be  between  0.3 Ad  and 
0.5 Ad-  Here  Ad  is  the  Debye  length, 


(3) 


where  €  is  the  permittivity,  Ab  is  the  Boltzmann  constant,  T 
is  the  temperature,  and  no  the  carrier  density. 

Stability  requires  that  the  fastest  likely  carrier  not  traverse 
a  full  grid  cell  over  the  duration  of  a  time  step,  At.  A  second, 
typically  more  restrictive  At  requirement  comes  from  the 
necessity  to  resolve  plasma  oscillations,  in  order  to  maintain 
charge  balance  across  the  testing  domain  and  avoid  charge 
instabilities  [1].  The  plasma  frequency  cov  is  defined  by 


(4) 


The  corresponding  time  step  requirement  is 
0.5 

At  < — .  (5) 

<z>p 

Note  that  these  stability  and  accuracy  restrictions  result  from 
the  incorporation  of  the  Poisson’s  equation  solver;  these  re¬ 
quirements  are  not  intrinsic  to  the  EMC.  If  one  is  interested 
in  resolving  transient-regime  transport,  At  should  addition¬ 
ally  be  much  smaller  than  the  relaxation  time  ( At  r). 

In  the  typical  EMC  implementation,  changes  in  carrier 
energy  and  momentum,  from  individual  scattering  events, 
are  calculated  dynamically  during  the  simulation.  These  cal¬ 
culations  are  simple  for  low-energy  electrons  where  the  en¬ 
ergy  is  parabolic  in  momentum,  as  well  as  for  intermediate- 
energy  electrons  where  the  simple  first-order  non-parabolic 
approximation  accurately  describes  the  band  structure  [2]. 
However,  these  simple  approximations  do  not  describe  high- 
energy  carrier  dynamics  well  [45].  Such  cases  require  the 
incorporation  of  the  full  band  structure  in  EMC,  known  as 
the  full-band  Monte  Carlo  [47-49].  This  technique  provides 
a  highly  accurate  description  of  high-field  transport  [44]  but 
the  added  computational  burden  is  significant. 

The  Cellular  Monte  Carlo  (CMC)  method,  originally  de¬ 
veloped  by  Kometer  et  at.  [50],  was  designed  to  reduce 
this  computational  burden,  but  the  CMC  method  itself  re¬ 
quires  a  sizable  amount  of  computer  storage  to  properly 
maintain  the  required  energy  and  momentum  conservation 
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laws.  A  hybrid  EMC/CMC  method  reduces  this  computa¬ 
tional  overhead,  and  offers  a  trade-off  in  terms  of  increased 
memory  usage.  The  particle-based  simulator  used  in  Sect.  4 
of  this  work  is  based  upon  the  hybrid  approach  developed 
by  Saraniti  and  Goodnick  [51].  The  general  algorithm  of 
the  full-band  device  simulator  follows  the  same  steps  as 
those  previously  shown  in  Fig.  1;  now  all  possible  initial 
and  final  states  for  every  scattering  mechanism  are  precom¬ 
puted  and  stored,  to  be  accessed  during  the  simulation.  Us¬ 
ing  the  hybrid  EMC/CMC  (henceforth  referred  to  as  CMC 
for  brevity),  bulk  material  simulations  have  been  demon¬ 
strated  to  run  50  times  faster  than  with  the  traditional  EMC 
method  [51]. 


2.2  Finite-Difference  Time-Domain  (FDTD)  simulation 


FDTD  is  a  highly  accurate  and  efficient  computational 
technique  for  modeling  electromagnetic  wave  interactions 
with  physical  structures.  Advances  in  absorbing  boundary 
conditions,  dispersive  and  nonlinear  materials  modeling, 
low-numerical-dispersion  schemes,  unconditionally  stable 
schemes,  and  incident  wave  source  conditions  make  FDTD 
a  highly  attractive  tool  for  electromagnetic  analysis.  Some 
details  on  these  advances  are  given  here,  and  further  infor¬ 
mation  is  available  in  Ref.  [5]. 

FDTD  is  a  direct  numerical  solution  of  the  time-depend¬ 
ent  Maxwell’s  curl  equations, 


d(fiH) 

dt 


=  -  V  x  E-M, 


d(€E) 

dt 


=V  xH-(j  +  crE) 


(6) 


where  E  and  H  are  the  electric  and  magnetic  fields,  respec¬ 
tively,  6,  /x  and  o  are  the  permittivity,  permeability  and  con¬ 
ductivity  of  the  medium,  respectively,  and  J  and  M  are  elec¬ 
tric  and  magnetic  source  current  densities.  Note  that  (6)  may 
be  expressed  as  six  coupled  scalar  equations,  each  a  partial 
differential  equation  involving  E ,  H ,  J ,  and  M  vector  com¬ 
ponents. 

The  fully  explicit  FDTD  algorithm  is  obtained  by  nu¬ 
merically  approximating  the  spatial  partial  derivatives  in  (6) 
with  centered  finite  differences,  and  numerically  integrat¬ 
ing  the  resulting  system  of  spatial  difference  equations  with 
respect  to  time  via  a  centered  finite-difference  approxima¬ 
tion  of  the  temporal  partial  derivatives.  The  staggering  of 
E  and  H  in  time  yields  an  efficient  leapfrog  time-marching 
scheme  wherein  all  of  the  E  components  are  updated  at  time 
step  n  (corresponding  to  a  physical  time  of  n  At)  using  pre¬ 
viously  stored  H  data,  and  then  all  of  the  H  components 
are  updated  at  time  step  n  +  1/2  using  the  just-computed 
E  data.  The  FDTD  method  does  not  require  any  formal 
treatment  of  Maxwell’s  divergence  equations  (e.g.  Gauss’s 


Fig.  2  (Color  online)  Illustration  of  the  staggered  electric  and  mag¬ 
netic  field  components  about  a  single  Yee  cell  in  a  3D  space  lattice. 
The  axis  origin  is  positioned  at  index  (i,  j,k ),  and  the  opposite  comer 
has  index  (i  +  1,  j  +  1,  k  +  1).  Each  electric  field  component,  indi¬ 
cated  in  red,  is  surrounded  by  four  circulating  magnetic  field  compo¬ 
nents.  Likewise,  each  magnetic  field  component,  indicated  in  blue,  is 
surrounded  by  four  circulating  electric  field  components 


Laws)  because  those  equations  are  implicitly  enforced  by 
Yee’s  space-time  grid,  as  shown  in  [5]. 

Figure  2  illustrates  the  staggered  sampling  of  E  and  H 
over  one  grid  cell  of  dimensions  Ax,  Ay,  and  Az,  in  a  3D 
Cartesian  spatial  lattice.  Each  vector  component  of  J  is  spa¬ 
tially  collocated  with  the  corresponding  component  of  E\ 
M  and  H  are  similarly  collocated.  Position  within  the  dis¬ 
crete  grid  is  defined  in  terms  of  grid  indices  (/,  j,k).  The 
grid  point  at  (/,  j,  k)  has  physical  location  [x (/ ) ,  y(j),  z(k)] 
where  x(i)  =  i  Ax,  y(j )  =  j  Ay,  and  z(k)  =  kAz. 

If  we  assume  all  structures  and  fields  to  be  z -invariant  so 
that  d/dz  ->  0,  we  obtain  a  simpler  two-dimensional  formu¬ 
lation.  The  2D  TEZ  FDTD  grid  describes  Hz,  Ex  and  Ey; 
for  most  media,  the  other  three  field  components  decouple 
from  these  three,  and  may  therefore  be  ignored  for  this  dis¬ 
cussion.  Note  that  the  TEZ  components  reside  on  the  x-y 
plane  in  Fig.  3. 
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Electromagnetic  energy  may  be  introduced  into  the 
FDTD  grid  through  the  application  of  nonzero  current  den¬ 
sities,  /  and  M,  which  may  be  assigned  any  temporal 
behavior  in  any  geometrical  configuration.  The  total-field 
scattered-field  (TFSF)  formulation  applies  electric  and  mag¬ 
netic  currents  along  a  transparent  closed  surface  within  the 
FDTD  grid  [5].  This  technique  calculates  J  and  M  ac¬ 
cording  to  surface  equivalence  [52]  to  produce  a  homoge¬ 
neous  propagating  plane  wave  within  the  closed  surface.  The 
analytic-field-propagation  (AFP)-TFSF  formulation  is  used 
in  Sect.  3  to  allow  plane  wave  interaction  with  conductive 
material  [53-55]. 

In  open-region  simulations,  grid  boundaries  are  truncated 
with  absorbing  boundary  conditions  such  as  the  perfectly- 
matched  layer  (PMF),  which  reduce  reflections  from  the 
grid  boundary  by  at  least  80  dB  [56,  57].  Thus,  the  main 
grid  is  electromagnetically  isolated  from  the  grid  boundary, 
so  that  a  finite  grid  may  model  an  infinite  space.  The  ma¬ 
jority  of  FDTD  implementations  described  in  Sects.  3  and  4 
use  convolutional  PMF  (CPMF)  [5,  58]. 

The  standard  FDTD  method  is  second-order-accurate  in 
both  space  and  time.  The  algorithm  exhibits  numerical  dis¬ 
persion  and  conditional  stability,  which  together  constrain 
the  choice  of  grid  cell  size  and  time  step.  Sub- wavelength 
sampling  in  space  is  required  to  properly  capture  the  high¬ 
est  spatial  frequencies  of  interest  in  the  problem  and  to  ade¬ 
quately  suppress  the  numerical  phase  velocity  errors.  Typi¬ 
cally,  10-20  samples  per  smallest  wavelength  of  interest  are 
needed  [5].  The  sampling  in  time  is  chosen  to  satisfy  the 
following  Courant  stability  condition  [59]: 

A  t< - .  1  (7) 

(x?)2  +  (^)2  +  (xi)2 

where  umax  represents  the  highest  wave  speed  in  the  sim¬ 
ulation  domain.  For  a  3D  lattice  where  Av  =  Ay  =  Az, 
At  <  Ax / (vmSLXV3) ,  whereas  for  a  2D  square  lattice,  At  < 
Av/(umaxV/2).  As  an  example,  consider  a  3D  simulation 
of  bulk  doped  silicon  at  room  temperature,  with  cr  =  11.8, 
Ax  =  20  nm,  and  n o  =  1016  cm-3.  The  Poisson  solver  time 
step  restriction  from  (5)  gives  A^p0iSS0n  <  0.29  ps.  At  the 
same  time,  FDTD  stability  requires  A^fdtd  <0.13  fs.  The 
maximum  Poisson  solver  time  step  is  more  than  three  orders 
of  magnitude  larger  than  the  maximum  FDTD  time  step. 

If  both  EMC  and  FDTD  are  updated  at  every  time  step, 
particle  positions  and  momenta  are  updated  far  more  of¬ 
ten  than  required,  at  the  cost  of  a  dramatic  increase  in  the 
computational  burden.  To  alleviate  the  burden  caused  by 
the  Courant  stability  limit  inherent  to  the  standard  FDTD 
technique,  we  use  the  alternating-direction-implicit  (ADI)- 
FDTD  technique  in  the  analyses  of  Sect.  4.  ADI-FDTD  is 
an  unconditionally  stable  algorithm;  that  is,  the  choice  of 
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Fig.  4  Flowchart  for  the  combined  EMC-FDTD  simulation  tool.  The 
carrier  transport  solver  (EMC  or  CMC)  acts  on  field  calculations  via 
the  spatially-varying  current  density  J.  The  electrodynamic  solver  acts 
on  charge  transport  calculations  via  spatially-varying  electric  and  mag¬ 
netic  fields  E,  H 

A t  is  no  longer  restricted  by  stability  requirements,  but  in¬ 
stead  by  much  looser  accuracy  requirements.  Details  on  this 
technique  are  available  in  Refs.  [5,  60-63]. 

2.3  Coupled  EMC-FDTD  simulation 

In  the  combined  EMC-FDTD  solver,  electric  and  magnetic 
fields  from  FDTD  influence  EMC  carrier  motion  accord¬ 
ing  to  the  Forentz  force,  (2).  Microscopic  currents  resulting 
from  carrier  motion  in  the  EMC  influence  FDTD-computed 
field  values  according  to  Maxwell’s  curl  equations,  (6). 
A  schematic  of  the  coupling  between  the  solvers  is  given 
in  Fig.  4. 

The  combination  of  a  full- wave  FDTD  solver  with  a  de¬ 
vice  simulator  based  upon  a  stochastic  transport  kernel  is 
conceptually  straightforward,  but  the  development  and  im¬ 
plementation  of  such  an  algorithm  can  be  challenging,  as 
we  will  see  with  specific  examples  in  the  following  sections. 
Simultaneous  adherence  to  the  stability  requirements  of  the 
Poisson  solver  and  FDTD  results  in  a  heavy  computational 
burden  in  the  combined  solver.  We  will  describe  two  differ¬ 
ent  approaches  to  this  issue. 

The  first  approach,  described  in  Sect.  3,  is  appropriate 
to  the  examination  of  bulk  materials  under  electromagnetic 
field  irradiation,  where  charge  distribution  is  expected  to 
be  spatially  constant.  If  doping  density  is  sufficiently  low 
(i no  <  1017  cm-3),  plasmon  interactions  impact  material  pa¬ 
rameters  only  negligibly,  and  the  Poisson’s  equation  solver 
is  not  necessary  for  the  calculation.  In  this  case,  the  stabil¬ 
ity  requirements  are  only  slightly  modified  from  those  of 
FDTD  [14],  and  the  computational  burden  is  substantially 
lightened.  This  approach  is  used  in  Sect.  3 

The  electrostatic  field  calculation  is  vital  to  accurate  de¬ 
vice  analysis.  In  microwave  devices  a  splitting  of  the  dc  and 
ac  components  and  a  simultaneous  solution  of  the  Poisson 
equation  for  the  dc  components  and  the  FDTD  solution  to 
the  ac  components  is  necessary.  Adherence  to  stability  crite¬ 
ria  of  both  solvers  prolongs  computation  times.  The  second 
approach,  described  in  Sect.  4,  relaxes  the  FDTD  stability 
requirements  on  A t  through  implementation  of  ADI-FDTD. 
Computational  burden  is  further  lifted  through  the  use  of  the 
Cellular  Monte  Carlo  (CMC)  method. 

The  EMC  tracks  carrier  positions  in  continuous  space, 
but  most  commonly  used  numerical  electromagnetics 
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solvers,  including  FDTD  and  Poisson’s  equation  solvers, 
calculate  electromagnetic  quantities  at  discrete  grid-based 
locations.  The  quasi-static  EMC-Poisson  solver  assumes  an 
instantaneously  electrostatic  system;  in  Maxwell’s  equations 
3/3 1  0.  If  we  neglect  the  effect  of  H  on  carriers,  the  only 

non-trivial  equality  is  Gauss’s  law, 

V-£  =  -,  (8) 

6 

where  p  is  charge  density,  which  is  incorporated  in  the  EMC 
through  the  Poisson  solver. 

The  EMC  and  the  electromagnetic  solver  must  interact 
through  a  continuous-to-discrete  spatial  mapping.  In  EMC- 
Poisson  simulations  this  is  called  the  charge  assignment 
scheme  (CAS)  [4,  64].  The  charge  density  associated  with 
each  carrier  in  continuous- space  is  assigned  to  discrete  grid 
locations  according  to  the  CAS.  The  Poisson’s  equation 
solver  uses  this  instantaneous  charge  density,  along  with  the 
known  permittivity,  to  find  the  grid-based  electrostatic  po¬ 
tential.  Partial  derivatives  of  the  potential  produce  the  grid- 
based  electric  field,  which  is  then  interpolated  in  continuous- 
space  to  determine  the  force  on  the  particle,  with  interpola¬ 
tion  weights  identical  to  those  in  the  CAS  [18,  64]. 

The  charge  assignment  scheme  is  the  backbone  of  the 
simulated  carrier-field  interaction.  Implementation  of  each 
scheme  may  be  more  or  less  appropriate  depending  on  lo¬ 
cal  simulation  particulars;  inappropriate  choice  of  CAS  may 
lead  to  non-physical  electric  fields  at  contacts  and  dielectric 
interfaces  [64] .  Inconsistent  coupling  may  cause  a  nonphys¬ 
ical  carrier  self-force  [4],  where  a  carrier  is  acted  on  by  its 
own  potential.  Minimization  of  these  effects  has  driven  ex¬ 
tensive  research  in  the  application  of  CAS  in  rectangular- 
grid  EMC-Poisson  solvers  [4,  64-66].  The  conceptual  ad¬ 
vances  resulting  from  this  work  have  been  invaluable  in 
establishing  the  Ensemble  Monte  Carlo  technique  as  the 
benchmark  against  which  other  computational  electronics 
methods  are  gauged. 

These  advances  were  made  with  the  assumption  of  a 
rectangular  grid,  with  field  components  positioned  at  grid 
cell  corners.  The  noncollocated  field  components  of  FDTD 
change  the  grid  geometry  sufficiently  that  it  cannot  be  taken 
for  granted  that  the  results  of  EMC-Poisson  research  apply. 
Instead  of  assigning  a  single  scalar  quantity,  charge,  to  a  sin¬ 
gle  grid,  we  must  now  assign  the  current  density  vector  com¬ 
ponents  to  each  of  the  three  noncollocated  grids.  Early  con¬ 
tributions  to  the  description  of  current  assignment  schemes 
are  given  in  Refs.  [9,  14,  18,  21]. 

The  electrodynamic  EMC-FDTD  solver  assumes  full 
time  dependence  of  the  fields,  and  incorporates  this  time 
dependence  through  Maxwell’s  curl  equations,  given  in  (6). 
This  does  not  mean  that  Gauss’s  law  no  longer  applies,  how¬ 
ever;  Gauss’s  law  is  indeed  satisfied  in  the  EMC-FDTD  by 


locally  enforcing  current  continuity,  given  by 


Current  continuity  is  easily  satisfied  in  the  typical  EMC- 
Poisson  implementation.  In  that  case,  p  and  J  are  calculated 
at  the  same  locations  on  the  grid,  and  as  a  result  current 
continuity  is  ensured  by  using  the  same  weighting  scheme 
for  the  two  quantities.  In  the  EMC-FDTD  implementation, 
however,  p  and  J  are  defined  at  separate  locations  on  the 
grid.  The  approximations  inherent  to  the  charge  and  current 
assignment  schemes  mean  that  current  continuity  must  be 
enforced  explicitly. 

Several  current  assignment  schemes  that  ensure  local 
charge  conservation  have  been  developed.  Marder  presented 
a  modified  version  of  Maxwell’s  equations  to  prevent  vio¬ 
lations  of  Gauss’s  law  in  numerical  tests  [67].  Villasenor 
and  Buneman’s  method  forced  exact  charge  conservation 
in  the  current  density  calculation  for  simply  shaped  parti¬ 
cles  [68].  A  generalized  charge-conserving  current  calcula¬ 
tion  technique,  for  particles  of  arbitrary  shape,  was  given  by 
Esirkepov  [69].  A  2005  comparison  showed  no  difference 
in  the  quality  of  results  between  the  Villasenor-Buneman 
method  or  the  Esirkepov  method  for  similar  particles  [70]. 
That  study  highlighted  the  necessity  of  local  charge  conser¬ 
vation  in  particle-in-cell  codes. 

Here,  we  describe  the  fundamental  EMC-FDTD  current 
assignment  scheme.  The  method  given  below  assumes  the 
Villasenor-Buneman  [68]  current  conservation  technique. 
The  reader  is  referred  to  Refs.  [68-70]  for  implementation 
details.  First  we  outline  the  EMC-Poisson  charge  assign¬ 
ment  schemes,  following  the  description  given  in  Ref.  [64], 
to  establish  notation  and  provide  a  basepoint  for  compari¬ 
son.  Then  we  detail  the  general  EMC-FDTD  current  assign¬ 
ment  scheme. 

2.3.1  EMC-Poisson  NGP  and  CIC  schemes 

The  two  most  commonly  used  charge  assignment  schemes 
are  the  nearest-grid-point  (NGP)  scheme  and  the  cloud-in¬ 
cell  (CIC)  scheme  [1,  65,  71].  The  NGP  scheme  treats  the 
carrier  as  a  8 -function  of  charge:  each  carrier’s  entire  charge 
is  assigned  to  the  nearest  grid  location. 

In  the  CIC  scheme,  each  carrier’s  charge  is  assumed  to  be 
evenly  distributed  over  a  finite  region  [65].  The  width  of  this 
charge  “cloud”  may  be  chosen  to  suit  the  implementation;  in 
this  discussion  we  will  assume  each  cloud  to  span  one  grid 
cell  in  each  direction.  The  CIC  scheme  assigns  each  carrier’s 
charge  density  among  the  grid  points  bounding  the  element 
containing  the  carrier;  see  Fig.  5.  The  various  charge  assign¬ 
ment  schemes  have  different  advantages,  and  each  is  appro¬ 
priate  under  different  conditions.  The  validity  and  accuracy 
of  these  schemes  and  others  were  examined  in  great  detail 
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Fig.  5  CIC  schematic.  The  dashed  lines  show  the  locations  of  grid 
lines,  and  the  dark  circles  show  grid  points.  The  open  circle  indicates 
a  carrier  located  at  (x,  y).  In  the  CIC  scheme,  the  carrier  is  assumed  to 
have  charge  spread  over  a  finite  region,  shown  by  the  grey  square.  In 
this  2D  example,  the  grid  element  containing  the  carrier  is  bounded  by 
four  grid  points,  with  indices  (i^,  j±) 


nearest  grid  point,  so  that  one  of  the  w(i,  j )  in  (10) 
has  value  one  and  the  rest  have  value  zero.  In  CIC,  we 
allow  i  =  floor(x/Ax)  and  j  =  floor (y/ Ay),  and  use 
;±  =  i  +  1/2  ±  1/2  and  J  =  j  +  1/2  ±  1/2  in  (10)  and 
(11).  In  both  schemes,  we  add  qw(i,  j)/V  to  pij  for  all 
(hj)- 

2.  Find  the  electrostatic  solution:  Starting  with  pij  and  6, 
Poisson’s  equation  gives  the  grid-based  potential  (pij. 

3.  Calculate  grid-based  E :  The  electric  field  is  calcu¬ 
lated  on  each  mesh  point  (/,  j)  as  Ef  j  =  —  ((pi+ij  — 

0i-l,;)/(2Ax)  and  Ejj  =  -(</>, j+ 1  -  <pij-i)/(2Ay). 

4.  Interpolate:  In  NGP,  the  field  acting  on  the  carrier  is 
assumed  to  be  that  at  the  nearest  grid  location.  In  CIC, 
E(x,  y)  is  found  using  the  w(i,  j )  calculated  for  step  1: 

Ex(x,y)=  J2  Mi,j)Ef  j 
i±J± 

and 


in  references  [64,  65].  Improvements  and  various  nuances  of 
implementation  were  described  in  Ref.  [64] . 

In  this  discussion,  we  examine  only  the  2D  case;  exten¬ 
sion  to  the  3D  mesh  is  straightforward.  Consider  a  carrier 
with  charge  q  positioned  at  (x ,  y )  in  continuous  space,  as 
shown  in  Fig.  5.  We  define  in  the  x-y  plane  a  corresponding 
uniform  mesh,  (/,  j),  where  i  e  [1,  A)]  and  j  e  [1,  Nj].  The 
grid  spacing  is  given  by  Ax  in  the  x -direction  and  Ay  in 
the  y-direction.  The  volume  of  one  grid  cell  is  V  =  l  Ax  Ay, 
where  l  is  the  assumed  domain  depth  in  the  z-direction. 

We  define  i+ ,  i~ ,  j+  and  j~  as  the  indices  of  the  grid 
points  at  the  comers  of  the  cell  that  contains  the  carrier 
(Fig.  5).  w(i,  j)  is  the  charge- assignment  weight  at  grid 
point  (/,  j),  and  the  grid  position  is  as  described  in  Sect.  2.2. 
For  the  CIC  scheme,  the  distribution  weights  are  given  by 

=  wxwy, 

W(/  +  ,  j~)  =  (1  -  WX)Wy, 

(10) 

W(l  ,j  +  )  =  WX(l~Wy), 

+  ,  J+)  =  (1  -  Wjc)(1  -  Wy), 

where 

x{i+)  —  x  v(/+)  —  v 
x(i+)-x(i~)  Ax  ’ 

y(7+)  -  y  y(i+)  -  y 

Wy  y(j+)-y(j~ )  Ay 

Note  that  charge  conservation  dictates  j±  w(i,  j)  =  1. 

1.  Charge  assignment:  The  point  charge  density  q/V  is 
assigned  to  the  grid,  producing  a  grid-based  charge  den¬ 
sity  pij.  In  NGP  the  entire  charge  q  is  assigned  to  the 


Ey(x,y)=  J2  w(iJ)Ejr 

The  NGP  scheme  is  computationally  simple,  but  inaccu¬ 
racies  are  introduced  because  the  NGP  scheme  effectively 
shifts  the  carrier  position  to  the  nearest  grid  location  for  all 
electromagnetic  interactions.  Additionally,  this  scheme  does 
little  to  counteract  local  charge  density  fluctuations  resulting 
from  a  finite  carrier  ensemble  size  [1,  65]. 

CIC  is  computationally  more  intensive  than  NGP,  but  it 
vastly  improves  smoothing  of  local  charge  density  fluctu¬ 
ations,  often  permitting  smaller  ensemble  sizes  and  a  de¬ 
creased  overall  computational  burden.  CIC  has  been  shown 
to  produce  nonphysical  fields  at  metal  contacts,  however, 
and  therefore  must  be  used  judiciously  [64,  65]. 

2.3.2  EMC-EDTD  NGP  and  CIC  schemes 


Consider  a  point  charge  with  the  current  density  contribu¬ 
tion  J  =  (Jx,  Jy)  located  at  (x,  y)  in  a  2D  FDTD  mesh,  as 
shown  in  Fig.  6.  The  computational  grid  shown  incorporates 
Hz,  Ex,  and  Ey  fields.  Because  the  Jx  and  Jy  components 
are  noncollocated  in  the  FDTD  formulation,  the  problem  is 
not  straightforward.  The  x-  and  y-current  densities  must  be 
considered  independently. 

Assume  a  CIC  scheme  with  a  rectangular  charge  cloud 
of  height  p  =  q/V.  A  carrier  moves  from  rn  =  (xw,  yn)  to 
rn+\  =  (xn+i,yw+i)  over  a  single  time  step  At.  For  now 
we  assume  rn  and  rn+\  are  in  the  same  grid  cell,  and  the 
cell  is  bound  by  indices  (/,  j)  and  (/  +  1,  j  +  1).  Then,  the 
Villasenor-Buneman  method  [68]  defines  the  current  density 
contribution  of  that  carrier  by 


r  Pi 

Ji+\J  =  M(Xn+1 


*n)  ^ 
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Fig.  6  The  carrier  at  (x ,  y)  in  the  staggered  FDTD  grid.  Only  the  Ex 
and  Ey  field  components  are  shown 
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(12) 


If  the  carrier  passes  through  a  grid  cell  boundary  (or  two 
boundaries,  if  the  carrier  is  located  at  the  corner  of  a  grid 
cell),  the  current  calculation  is  split  into  two  (or  three)  parts, 
where  each  part  is  used  to  find  the  carrier’s  current  contribu¬ 
tion  within  each  cell.  See  Refs.  [68,  70]  for  further  details. 
If  J  is  defined  in  this  manner,  current  continuity  is  ensured. 


1.  Current  assignment:  The  current  density  of  the  parti¬ 
cle  is  calculated  from  the  change  in  the  particle  position 

over  At.  J*+1/2j>  1/2, j+l’  Jlj+ 1/2  and  Ji+\,j+\/2  are 

calculated  according  to  (12)  and  the  procedures  given  in 
Refs.  [68-70]. 

2.  Find  the  electrodynamic  solution:  Using  Jx ,  Jy  and  e, 
FDTD  gives  the  electric  and  magnetic  fields  Ex ,  Ey  and 
Hz. 

3.  Interpolate  field:  It  is  necessary  to  spatially  average 
the  electric  fields  so  that  they  are  known  at  the  grid 
intersections,  in  order  to  avoid  particle  self-force  [70]. 

Kj  =  (E+H2J  +  £f-l/2j)/2  and  Ejj  =  (EfJ+1/2  + 
Efj_ i/2)/2.  The  electric  field  at  the  carrier’s  location 
(x,  y)  is  found  by  interpolating  between  the  nearest  grid 
points  following  the  same  procedure  as  that  used  in  the 
CIC  charge  assignment  scheme. 


3  THz  conductivity  of  doped  bulk  silicon 

The  plasma  frequency  and  characteristic  scattering  rate  of 
moderately  doped  semiconductors  typically  fall  within  the 
THz  regime  [72].  For  this  reason,  semiconductor  carrier  re¬ 
sponse  to  THz-frequency  stimulation  depends  strongly  on 
the  specific  form  of  the  material  transport  parameters.  The 
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bulk  materials  behavior  in  the  THz-range  is  extremely  sen¬ 
sitive  to  small  deviations  in  doping  density  and  tempera¬ 
ture  [73,  74]. 

In  this  section,  we  explore  the  use  of  the  global  EMC- 
FDTD  solver  in  the  context  of  doped  bulk  silicon  materials 
characterization.  First,  for  reference,  we  review  THz-regime 
experimental  characterizations  of  doped  silicon  in  Sect.  3.1. 
Then  in  Sect.  3.2  we  describe  the  specific  coupled  EMC- 
FDTD  solver  used  to  simulate  THz-frequency  electromag¬ 
netic  plane  wave  interactions  with  doped  bulk  silicon,  and 
the  method  we  use  to  extract  the  frequency-dependent  ef¬ 
fective  conductivity  from  the  computed  fields  and  currents. 
We  examine  simulation  performance  in  the  context  of  pre¬ 
dicted  conductivity  convergence  under  variation  in  grid  cell 
size,  ensemble  size,  averaging  technique,  and  impedance 
mismatch.  Finally,  we  compare  the  simulation  results — that 
is,  the  predicted  effective  conductivity — with  published  ex¬ 
perimental  results  for  doped  silicon  at  THz  frequencies. 

3.1  Experimental  characterization 

Since  the  advent  of  THz  time-domain  spectroscopy  (THz- 
TDS)  20  years  ago,  extensive  experimental  work  has  exam¬ 
ined  the  THz  characteristics  of  doped  bulk  silicon  [72-84]. 
Pure  undoped  silicon  is  almost  completely  transparent  and 
nondispersive  under  THz-frequency  radiation,  more  so  than 
quartz,  sapphire,  or  fused  silica,  making  silicon  a  very  attrac¬ 
tive  THz-optics  material  [73].  However,  this  behavior  is  only 
observed  in  extremely  pure  samples.  Carriers  introduced  by 
low-level  doping  and  impurities  strongly  impact  the  THz- 
frequency  materials  characteristics  of  silicon  [73].  As  dop¬ 
ing  density  is  increased,  silicon  becomes  quite  opaque  to 
THz  radiation  [72].  Because  pure  silicon  is  extremely  trans¬ 
parent,  it  is  reasonable  to  assume  that  carrier-field  dynamics 
within  doped  silicon  are  solely  responsible  for  doped  sili¬ 
con’s  opacity  under  THz  radiation  [73]. 

The  seminal  THz-TDS  characterizations  extracted  the 
complex  transmission  spectrum  of  semiconductor  samples 
through  the  use  of  freely -propagating  THz  beams.  The  high 
accuracy  of  this  technique  is  mainly  limited  by  uncertainties 
in  sample  thickness  that  lead  to  minor  uncertainties  in  the 
transmission  amplitude,  and  comparatively  large  (~10%) 
uncertainties  in  the  phase  [72]. 

As  carrier  density  increases  to  around  10 17  cm-3,  the 
increased  opacity  of  silicon  precludes  transmission-based 
characterization  [79].  The  development  of  reflection-based 
THz-TDS  in  1996  [85]  has  since  permitted  extensive  study 
of  optically  dense  media  [74,  79].  In  reflecting  THz-TDS 
the  exact  placement  and  orientation  of  the  sample  is  crit¬ 
ically  important.  In  analogy  to  uncertainties  from  sample 
thickness  in  transmission  THz-TDS,  here  small  deviations 
in  sample  position  contribute  to  uncertainty  in  the  reflection 
spectrum  amplitude  and  larger  uncertainties  in  phase.  This 
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Table  1  Summary  of  experimental  results  on  doped  bulk  silicon  under 
THz  radiation,  as  obtained  from  THz-TDS,  chronologically  ordered. 
Values  in  bold  were  calculated  for  comparison  purposes,  based  on  data 
from  the  referenced  work.  The  table  shows  the  conductivity  model  as¬ 
sumed  in  the  original  work;  the  material  dc  characteristics  including 

resistivity,  doping  density,  and  mobility;  the  fitted  THz  characteristics 
including  doping  density  and  mobility;  the  percent  change  from  dc  to 
THz  doping  density  and  mobility  values,  and  the  data  citation.  In  most 
cases,  both  no  and  /x  are  taken  as  fitting  parameters;  Ref.  [82]  assumes 
known  /x  for  all  calculations. 

Model 

dc 

THz 

Ano  % 

A/x  % 

Ref. 

p  (£2  cm) 

n0  (cm  3) 

/x  (cm2/V  s) 

n0  (cm-3) 

/x  (cm2/V  s) 

Drude 

1.15 

4.6  x  1015 

1302 

3.3  x  1015 

1680 

28 

-29 

[77] 

Drude 

8.1 

7.0  x  1014 

1349 

4.2  x  1014 

1820 

40 

-35 

[72] 

Drude 

8.15 

5.2  x  1014 

1352 

3.4  x  1014 

2340 

34 

-73 

[78] 

Drude 

0.21 

3.2  x  1016 

1083 

2.1  x  1016 

1420 

34 

-31 

[78] 

Cole-Davidson 

8.15 

5.2  x  1014 

1352 

4.0  x  1014 

2000 

23 

-47 

[78] 

Cole-Davidson 

0.21 

3.2  x  1016 

1083 

2.3  x  1016 

1280 

28 

-18 

[78] 

Drude 

1.1 

3.8  x  1015 

1350 

4.9  x  1015 

1350 

-29 

0 

[82] 

Drude 

2.5 

1.7  x  1015 

1350 

1.8  x  1015 

1350 

-6 

0 

[82] 

Drude 

10.0 

4.2  x  1014 

1350 

6.7  x  1014 

1350 

-59 

0 

[82] 

Drude 

1.1 

4.1  x  1015 

1306 

3.7  x  1015 

1560 

12 

-19 

[74] 

Drude 

0.14 

4.0  x  1016 

1036 

3.5  x  1016 

1280 

14 

-17 

[80] 

uncertainty  may  be  minimized  by  analytically  sweeping  the 
relative  positions  of  sample  and  reference  plane;  see,  for  ex¬ 
ample,  Ref.  [79].  While  careful  implementation  of  reflect¬ 
ing  THz-TDS  produces  excellent  results,  the  sensitivity  of 
the  technique  has  prompted  research  into  other  THz-regime 
characterization  methods  [81-84]. 

In  summarizing  available  experimental  results  on  the 
THz  characteristics  of  doped  silicon,  we  consider  only  those 
reports  which  provide  sufficient  data  for  comparison:  the 
room  temperature  dc  resistivity  or  assumed  equivalent  dop¬ 
ing  density,  as  well  as  THz-range  fitting  doping  density  or 
plasma  frequency  [72,  74,  77-81].  These  data  are  shown  in 
Table  1.  In  nearly  every  case,  these  results  were  extracted 
from  observed  reflection  and  transmission  spectra  assuming 
Drude  model  behavior,  where  the  Drude  model  conductivity 
is  given  by 


r  is  the  characteristic  carrier  scattering  time  and  cov  is  the 
plasma  frequency.  If  we  consider  only  the  data  in  Table  1 
in  which  carrier  density  and  mobility  are  taken  as  fitting 
parameters,  the  fit  to  the  Drude  model  requires  assumed 
doping  densities  that  average  ~27%  lower  than  those  cal¬ 
culated  from  the  known  dc  resistivity,  and  calculated  mo¬ 
bilities  ~32%  higher  than  those  found  at  dc.  This  sug¬ 
gests  that  r  energy-dependence  sufficiently  impacts  THz- 
frequency  properties  that  it  must  be  accounted  for  by  the 
conductivity  model  [72].  Several  other  models  have  been  de¬ 
veloped  to  incorporate  t  energy  dependence,  with  promising 
results  [72,  78]. 


The  ill-fit  of  the  Drude  model  to  doped  silicon  at  THz  fre¬ 
quencies,  and  the  absence  of  another  well-accepted  model, 
makes  prediction  of  THz-frequency  materials  properties  of 
arbitrarily  doped  silicon  difficult.  As  demonstrated  below, 
the  EMC-FDTD  method  has  the  potential  to  provide  detailed 
materials  characterization  of  doped  silicon  at  THz  frequen¬ 
cies. 

3.2  Simulation  domain 

For  both  EMC  and  FDTD  we  use  two-dimensional  (2D) 
computational  domains  defined  in  the  v-y  plane  (Fig.  8). 
The  EMC  domain  is  filled  with  (001)  doped  silicon  with 
doping  density  no  =  1017  cm-3.  The  typical  simulation  car¬ 
rier  ensemble  size  is  O(105).  To  permit  examination  of 
the  interaction  between  these  carriers  and  propagating  THz- 
frequency  electromagnetic  plane  waves,  we  embed  the  EMC 
domain  into  an  FDTD  domain. 

The  FDTD  simulation  testbed  is  a  semi-infinite  half  space 
of  doped  silicon  (Regions  B  and  C  in  Fig.  8)  and  a  semi¬ 
infinite  half  space  of  air  (Region  A  in  Fig.  8).  Region  A 
is  assigned  a  dielectric  constant  of  1  and  zero  conductiv¬ 
ity.  A  dielectric  constant  of  11.8  is  assigned  to  Regions  B 
and  C.  Region  C  is  filled  with  an  assumed  value  for  the  con¬ 
tinuous  bulk  dc  conductivity  of  doped  silicon,  d .  Region  B 
contains  the  embedded  EMC  domain;  the  current  J  in  (6)  is 
generated  solely  from  the  EMC. 

The  EMC  formulation  accounts  for  the  smooth  material 
interface  by  enforcing  specular  reflection  of  carriers  at  the 
left  and  right  boundaries  of  the  domain.  The  top  and  bot¬ 
tom  boundaries  are  given  periodic  boundary  conditions,  to 
allow  unrestricted  carrier  motion  in  the  vertical  direction. 
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Fig.  7  Flowchart  for  the  combined  EMC-FDTD  simulation  modeling 
tool 

Care  must  be  taken  at  the  edges  of  the  coupled  region  to  en¬ 
sure  that  EMC  current  density  does  not  extend  beyond  the 
coupled  region  for  any  charge  assignment  scheme. 

To  allow  finite-grid  representation  of  an  infinite  space, 
we  treat  FDTD  domain  boundaries  with  CPML  absorbing 
boundary  conditions  [5,  58].  Outward-propagating  waves 
that  reflect  from  the  domain  boundary  are  attenuated  within 
the  CMPL  by  more  than  80  dB,  so  that  the  main  grid  is 
electromagnetically  isolated  from  the  domain  boundary  to 
very  good  approximation.  We  assume  a  TEZ  mode  for  the 
electromagnetic  wave.  Plane  waves  are  introduced  with  a 
propagation  direction  normal  to  the  interface  via  the  AFP- 
TFSF  formulation  [53-55].  In  this  polarization,  the  incident 
plane  wave  will  excite  y -directed  currents  within  the  EMC- 
coupled  region.  The  periodic  boundaries  on  the  top  and  bot¬ 
tom  edges  of  the  EMC  domain  and  the  uniform  incident 
field  allow  bulk  EMC  simulation  within  the  finite  domain. 
Figure  7  shows  the  full  simulation  flowchart  for  a  coupled 
EMC-FDTD  solver. 

Figure  8  shows  the  2D  spatial  distribution  of  the  Ey 
phasor  amplitude,  extracted  via  discrete  Fourier  transform 
over  several  periods  of  electromagnetic  wave  oscillation. 
The  field  amplitude  is  attenuated  as  a  function  of  depth  in 
both  Regions  B  and  C.  Amplitude  decay  in  region  C  corre¬ 
sponds  to  that  expected  for  a  dielectric  with  conductivity  b . 
No  conductivity  is  enforced  in  the  coupled  Region  B ;  field 
decay  in  this  region  results  from  the  carrier-field  interaction, 
exhibiting  the  macroscopic  phenomenon  of  the  skin  effect. 

Noise  in  the  Ey  phasor  amplitude  in  Region  B  is  caused 
by  thermal  electron  motion;  Ex ,  Jx  and  Jy  also  exhibit  mi¬ 
nor  fluctuations  in  phase  and  amplitude.  Transient  fields, 
continually  sourced  by  this  noise,  propagate  to  the  grid 
boundary  and  are  attenuated  in  the  CPMF.  Without  high- 
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Fig.  8  Amplitude  of  Ey  field  phasor  extracted  by  DFT  over  several 
periods  of  electromagnetic  wave  oscillation,  where  white  corresponds 
to  high  field  amplitude  and  black  corresponds  to  low  field  amplitude. 
The  materials  interface  is  indicated  by  the  vertical  gray  line.  The  sec¬ 
tion  of  the  domain  shown  here  lies  within  the  AFP-TFSF  boundary; 
an  incident  plane  wave  is  sourced  from  the  left  boundary.  Region  B  is 
the  EMC-FDTD  coupled  region,  whose  boundary  is  marked  by  a  solid 
gray  box,  and  Regions  A  and  C  are  pure-FDTD  doped  silicon 


quality  absorbing  boundary  conditions  these  simulations 
would  suffer  from  the  same  continually  increasing  noise 
level  that  affected  the  early  simulations  described  in  Sect.  1 . 
To  reduce  the  impact  of  this  noise  on  conductivity  calcula¬ 
tions,  we  take  spatial  averages  over  small  regions  surround¬ 
ing  each  grid  location  in  the  extracted  phasor  quantities.  In¬ 
creasing  the  size  of  the  averaging  regions  decreases  pha¬ 
sor  quantity  noise.  The  averaging  region  size  must  not  be 
increased  beyond  ~l/20th  of  the  smallest  electromagnetic 
feature  of  interest. 

These  noise-reduced  phasor  quantities  are  used  in  the  ef¬ 
fective  linear-regime  conductivity  calculation, 


or(co)  = 


E(co)  •  J*(co) 
\E(co)\2 


(14) 


The  real  part  of  a  corresponds  to  power  dissipation  and 
is  evidenced  as  the  phasor  amplitude  decay  in  Fig.  8.  The 
imaginary  part  of  b  corresponds  to  phase  shift  between  E 
and  J  resulting  from  the  delay  in  material  response  to  ap¬ 
plied  fields.  As  the  electromagnetic  oscillation  frequency 
approaches  the  carrier  scattering  rate,  we  expect  the  delay 
in  material  response  to  applied  fields  to  increase. 

Figure  9  shows  the  variance  of  \Ey  \  in  Region  B  as  a 
function  of  electron  ensemble  size  Ne  for  several  grid  cell 
sizes  Ax.  As  expected,  larger  carrier  ensembles  show  de¬ 
creased  phasor  quantity  noise,  at  the  cost  of  an  increased 
computational  burden.  Figure  9  shows  dramatic  noise  reduc¬ 
tion  for  increased  Av .  This  improvement  directly  contrasts 
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Fig.  9  Average  noise  of  Ey  in  the  EMC-FDTD  coupled  region  for 
several  values  of  electron  ensemble  size  and  grid  cell  size.  E^01se  is 

calculated  as  <J(\Ey\2)  —  \{Ey)\2  and  is  normalized  by  \Ey\  at  the  ma¬ 
terial  interface.  Axref  ~  400  nm  is  a  reference  grid  cell  size.  Increasing 
the  size  of  the  ensemble  reduces  noise,  as  expected,  but  increasing  the 
grid  cell  size  produces  a  much  larger  improvement 


Fig.  10  Extracted  conductivity  with  varied  averaging  region  size,  elec¬ 
tron  ensemble  size,  and  surrounding  bulk  conductivity.  Nq  is  a  ref¬ 
erence  ensemble  size,  typically  O(105).  Larger  averaging  region  size 
leads  to  convergence  in  b .  Level  of  impedance  mismatch  between  Re¬ 
gions  B  and  C,  examined  through  modified  <r,  is  associated  with  line 
type.  See  text  for  discussion 


the  EMC-Poisson  solver  accuracy  requirements,  which  fa¬ 
vor  smaller  grid  cells  for  improved  electrostatic  calculations. 
The  EMC-FDTD  solver  achieves  significant  noise  reduction 
for  larger  Ax  by  including  a  larger  number  of  carriers  in 
each  Jx  and  Jy  grid  point  calculation. 

Figure  10  shows  b  as  a  function  of  phasor  quantity  av¬ 
eraging  region  size,  for  several  values  of  ensemble  size  and 
several  levels  of  impedance  mismatch  between  Regions  B 
and  C.  As  the  averaging  regions  size  is  increased,  b  con¬ 
verges.  Increased  ensemble  size  also  leads  to  convergence 
in  <7,  indicating  correspondence  between  decreased  phasor 
quantity  noise  and  convergence  in  b . 

Finally,  Fig.  10  allows  examination  of  the  impact  of  im¬ 
pedance  mismatch  between  the  EMC-FDTD  coupled  Re¬ 
gion  B  and  the  surrounding  pure-FDTD  Region  C.  In  any 
single  test,  the  conductivity  b  defined  in  Region  C  may  not 


match  the  conductivity  exhibited  by  the  EMC-FDTD  cou¬ 
pled  region.  The  resulting  impedance  mismatch  will  cause 
waveguiding  and  back-reflections  within  the  coupled  region. 
We  tested  the  effect  of  these  reflections  on  b  by  comparing 
extracted  conductivity  values  for  several  tests  where  b  was 
varied  through  three  values,  spanning  b .  Figure  10  shows 
the  results  of  this  experiment;  b  is  insensitive  to  impedance 
mismatch  between  Regions  B  and  C. 

3.3  Comparison  with  experiment 

We  compare  the  complex  b  with  experimental  results  ob¬ 
tained  via  reflecting  THz-TDS  [78].  The  dc  resistivity  of  the 
n-type  doped  silicon  sample  used  is  given  as  8. 15  Q  cm,  cor¬ 
responding  to  no  ~  5.5  x  1014  cm-3  [86].  This  doping  den¬ 
sity  is  used  in  the  EMC-FDTD  solver,  and  the  pure-FDTD 
a  is  modified  to  match  the  corresponding  dc  conductivity. 

The  analytical  best  fit  of  the  Cole-Davidson  model  to  ex¬ 
perimental  data  is  given  by  Jeon  and  Grischkowsky  [78]. 
The  Cole-Davidson  fit  can  be  regarded  as  a  faithful  rep¬ 
resentation  of  the  experimental  data.  Figure  11  compares 
the  EMC-FDTD-extracted  doped-silicon  conductivity  to  the 
best  fit  to  experimental  results.  The  global  solver  provides 
results  which  agree  well  with  experiment.  The  real  part  of 
the  conductivity  shows  excellent  agreement,  especially  at 
higher  frequency.  The  imaginary  component  of  the  conduc¬ 
tivity  shows  good  agreement,  even  though  small  phase  errors 
in  the  computationally  and  experimentally  extracted  con¬ 
ductivities  would  show  up  here.  In  the  numerical  tests,  no 
is  sufficiently  low  that  the  binary  Coulomb  interactions  be¬ 
tween  individual  carriers  should  not  substantially  contribute 
to  the  observed  conductivity,  but  because  the  plasma  fre¬ 
quency  cov  is  O(THz),  plasmon  scattering  may  be  important. 

This  preliminary  work  on  a  combined  EMC-FDTD 
solver  demonstrates  its  utility  in  characterization  of  doped 
semiconductors  at  THz  and  sub-THz  frequencies.  We  have 
shown  that  this  global  numerical  technique  provides  conduc¬ 
tivity  predictions  that  match  well  to  experiment.  The  EMC- 
FDTD  solver  exhibits  improved  accuracy  for  larger  Ax ,  in 
contrast  to  traditional  EMC-Poisson  accuracy  requirements. 
The  EMC-FDTD  global  solver  shows  promise  as  a  method 
for  full  characterization  of  doped  silicon  at  THz  frequencies. 


4  Global  modeling  of  microwave  transistors 

Microwave  semiconductor  device  analysis  in  the  small- 
signal  regime  requires  a  full  description  of  the  quasi-static 
dc  field  distribution  prior  to  the  application  of  a  small 
ac  modulation.  In  the  global  model  developed  recently  by 
Ayubi-Moak  et  al.  [37-41],  dc  device  behavior  is  obtained 
first,  as  the  steady-state  result  of  the  CMC-Poisson  simula¬ 
tion  at  a  given  bias  point.  The  calculated  dc  current,  /do  is 
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Fig.  11  THz  conductivity  of  doped  bulk  silicon.  The  EMC-FDTD 
conductivity  ( circles )  agrees  well  with  the  experimental  data  of 
Ref.  [78]  (presented  here  via  a  faithful  Cole-Davidson  analytical  fit, 
solid  curve ) 


used  to  source  the  dc  Maxwell’s  curl  equations, 

V  x  Edc  =  0,  (15a) 

Vxi  =  /dc.  (15b) 


An  ac  excitation  source  is  then  applied  to  the  device  via 
a  source  plane  approach  described  in  the  next  section.  The 
resulting  time- varying  distributions  of  the  EM  fields  are  then 
computed  by  solving  the  discretized  Maxwell’s  equations 
over  the  entire  FDTD  grid  (repeated  here  for  convenience) 
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Fig.  12  Flowchart  of  coupled  CMC/FDTD  device  simulator  [38] 


where  /t0 tal  is  now  the  total  current  density  (i.e.  the  ac  part 
plus  the  dc  part)  computed  by  the  CMC  solver  after  a  spec¬ 
ified  number  of  particle  free-flight  time  steps.  In  (16a)  and 
(16b),  a  lossless  medium  has  also  been  assumed  for  simplic¬ 
ity.  Note  that  during  the  first  FDTD  time  step  of  the  full- 
wave  simulation  in  which  the  ac  excitation  is  applied,  /t0 tal 
is  exactly  equal  to  /dc. 

Following  the  individual  update  of  each  field  component, 
the  new  ac  fields  are  then  added  to  the  stored  dc  fields  and 
the  updated  total  field  is  passed  back  into  the  CMC  solver 
and  used  to  compute  the  total  Forentz  force,  evolve  the  par¬ 
ticle  velocities,  and  recalculate  the  total  current.  Once  com¬ 
puted,  the  new  current  density  in  each  grid  cell  is  passed 
back  to  the  FDTD  solver  at  the  next  full-wave  time-step 
and  the  process  is  repeated  for  the  remainder  of  the  full- 
band/full-wave  simulation,  as  illustrated  in  Fig.  12. 

A  critical  component  in  the  full-wave  simulation  is  the 
manner  in  which  small- signal  voltage  perturbations  are  ap¬ 
plied  to  the  electrodes  on  the  FDTD  grid.  This  can  be  accom¬ 
plished  by  applying  the  proper  transverse  electric  field  dis¬ 
tribution  to  a  source  plane  at  the  front  end  of  the  simulation 
domain.  The  applied  electric  field  distributions  correspond 
to  the  solution  of  Poisson’s  equation  over  a  2D  source  plane 
for  a  set  of  applied  contact  voltages.  This  static  field  distrib¬ 
ution  is  then  modulated  by  the  appropriate  time  signature  of 
the  desired  excitation  function.  The  resulting  distributions 
are  then  applied  as  soft  sources  to  the  corresponding  field 
components  at  each  FDTD  time  step  during  the  full-wave 
portion  of  the  total  simulation. 

The  source  plane  is  located  along  a  cross-section  of  the 
3D  structure,  as  illustrated  in  Fig.  13,  for  a  simple  coplanar 
waveguide  (CPW)  geometry.  The  FDTD  grid  is  extended 
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Fig.  13  Example  of  simulation  domain  used  in  a  3D  CMC/FDTD  sim¬ 
ulation  showing  passive  and  active  sections 

both  in  the  front  of  and  to  the  back  of  the  actual  active  de¬ 
vice  region.  Both  the  passive  front-end  and  back-end  sec¬ 
tions  function  as  lossless  transmission  lines.  The  passive 
front-end  section  allows  the  sourced  excitation  wave  to  de¬ 
velop  into  the  proper  waveguide  mode  prior  to  its  interac¬ 
tion  with  the  active  device.  Both  passive  sections  are  defined 
as  extensions  of  the  intrinsic  undoped  semiconductor  mater¬ 
ial  beyond  the  physically  modeled  electrode  terminations  in 
the  CMC  simulator.  The  passive  sections  separate  the  active 
CMC  device  model  from  the  CPML  layers  (not  shown  in 
figure)  which  truncate  the  domain.  The  domain  region  above 
the  active  device  and  passive  regions  is  filled  with  free  space 
and  separates  the  top  of  the  circuit  from  the  upper  absorbing 
boundary  layer.  The  excitation  source  in  Fig.  13  is  applied  to 
the  gate  electrode  of  a  simple  three-terminal  device.  In  this 
illustration,  the  source  electrode  is  grounded  and  the  drain 
is  left  floating  (i.e.  not  grounded).  The  excitation  signal  is 
shown  being  applied  to  the  gate.  Not  shown  in  the  figure  is 
the  ground  plane  located  directly  below  the  substrate. 

4.1  Ultrafast  device  simulation  and  results 

The  first  direct  measurement  of  femtosecond  velocity  over¬ 
shoot  in  GaAs  and  InP  was  performed  by  Leitenstorfer  et 
al.  [87],  by  capturing  transient  THz  radiation  produced  by 
optically-generated  electron-hole  pairs  in  pin  diodes.  This 
work  motivated  the  subsequent  numerical  results  of  the  tran¬ 
sient  response  in  GaAs  and  InP  by  Wigger  et  al.  [88],  which 
used  CMC/Poisson  transient  simulation  (basically,  the  tran¬ 
sient  electric  field  was  assumed  to  satisfy  Poisson’s  equation 
at  every  point  in  time,  and  magnetic  field  was  neglected). 
Here  we  present  the  application  of  the  coupled  CMC/FDTD 
to  the  simulation  and  capture  of  radiated  EM  fields  emanat¬ 
ing  from  ultrafast,  high-frequency  devices  [37]. 

The  test  structure  was  a  3D  slice  of  intrinsic  GaAs 
500  nm  long  and  100  nm  wide  in  both  transverse  dimen¬ 
sions.  Contacts  were  simulated  on  each  end  of  the  structure 
allowing  for  the  necessary  reverse  bias  across  the  pin  diodes 
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Fig.  14  Cross-section  of  the  structure  used  in  the  THz  simulations 

used  in  the  original  experiments  in  [87].  A  500  nm  long  air 
region  was  added  to  the  left  end  of  the  structure  and  the  en¬ 
tire  simulated  domain  was  surrounded  by  a  split-field  per¬ 
fectly  matched  layer  (PML)  absorbing  boundary  region  ten 
cells  thick  in  each  direction.  A  500  nm  long  square  filament 
with  a  cross-sectional  area  of  400  nm2  was  centered  inside 
the  GaAs  region.  This  filament  represented  the  region  over 
which  simulated  electron-hole  pairs  were  injected  into  the 
device  structure.  Figure  14  shows  a  2D  cross-section  of  the 
test  structure.  This  device  was  simulated  using  a  uniform 
mesh  of  10  nm  spaced  grid  cells.  This  spatial  dimension  al¬ 
lowed  for  a  maximum  time  step  of  only  0.0016  fs  and  was 
the  limiting  factor  in  these  simulations. 

Two  sets  of  simulations  were  performed  on  the  test  struc¬ 
ture.  Both  simulations  were  run  for  a  total  of  350  fs.  In  each 
case,  25,000  electron-hole  pairs  were  injected  into  the  fila¬ 
ment  region  of  the  device  with  a  peak  injection  concentra¬ 
tion  of  5  x  10“ 14  cm-3  occurring  20  fs  after  the  start  of 
the  simulation.  In  the  first  experiment,  a  bias  of  —5  V  was 
applied  across  the  device,  corresponding  to  a  uniform  elec¬ 
tric  field  of  100  kV/cm.  In  the  second  experiment,  the  bias 
voltage  was  decreased  to  —1.2  V,  corresponding  to  a  uni¬ 
form  electric  field  of  24  kV/cm.  Because  each  simulation 
demonstrated  a  similar  set  of  results,  only  the  details  of  the 
100  kV/cm  simulation  are  presented  and  discussed  here. 

During  runtime,  snapshots  of  the  electric  and  magnetic 
field  components  were  taken  over  a  cross  section  of  the  y-z 
plane  at  distances  of  50  nm  (Fig.  15)  and  480  nm  (Fig.  16) 
from  the  left  contact.  These  points  were  chosen  in  order 
to  capture  the  near  and  far-field  characteristics  of  the  elec¬ 
tric  fields  surrounding  the  device.  In  Fig.  15,  a  snapshot 
of  the  Ey  field  was  taken  approximately  1.25  fs  after  the 
peak  of  the  injection  pulse.  It  is  clear  from  the  contour  plot 
that  the  field  is  dipolar  in  nature.  This  behavior  is  actually 
the  expected  response  from  an  elemental  oscillating  elec¬ 
tric  or  Hertzian  dipole  in  which  the  magnitude  of  the  elec¬ 
tric  field  intensity  is  directly  proportional  to  the  dipole  mo¬ 
ment  in  the  longitudinal  direction.  Under  the  influence  of  a 
high,  uniform  electric  field,  electrons  and  holes  are  acceler¬ 
ated  in  opposite  directions,  creating  an  electric  dipole  with 
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Fig.  15  Snapshot  of  Ey  taken  50  nm  from  the  left  contact,  1.25  fs 
after  the  peak  of  the  injection  pulse.  The  observed  dipolar  nature  of  the 
pulse  is  characteristic  of  the  near-field  region 
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Fig.  16  Snapshot  taken  480  nm  from  left  contact  of  transverse  Ey 
component  of  the  electric  field  74  fs  after  the  peak  of  the  injection 
pulse.  This  wavefront  resembles  a  plane  wave  (TEM  wave)  in  the  far 
field  region 


a  constantly  increasing  dipole  moment.  This  displacement 
of  both  electrons  and  holes  versus  simulation  time  is  plotted 
in  Fig.  17  and  validates  the  notion  of  an  oscillating  dipole. 
When  one  moves  into  the  far-held  region  (Fig.  16)  the  wave- 
front  becomes  TEM  or  plane  wave-like  in  nature.  The  re¬ 
sulting  drift  velocities  of  both  electrons  and  holes  are  shown 
in  Fig.  18.  Peak  overshoot  velocities  of  7.4  x  107  cm/s  for 
electrons  and  —1.5  x  107  cm/s  for  holes  were  determined 
and  are  in  good  agreement  with  those  reported  in  [88]  using 
just  the  CMC/Poisson  device  simulator. 


Time  [s] 


Fig.  17  Displacement  of  electrons  and  holes  in  the  device  versus  time 
for  the  simulation  run  at  100  kV/cm 


Fig.  18  Drift  velocities  of  electrons  and  holes  versus  time  for  the  sim¬ 
ulation  run  at  100  kV/cm 

4.2  Single  0.1  pm  gate  GaAs  MESFET 

The  coupled  device  simulator  has  also  been  applied  to  the 
simulation  of  a  simple  3D  MESFET.  In  [38],  a  0.1  pm, 
single-gate,  3D  MESFET  was  modeled  using  a  version  of 
the  coupled  simulator  which  used  a  conventional  3D  FDTD 
with  split-held  PML  absorbing  boundary  conditions.  Fourier 
decomposition  was  used  to  study  the  small  signal  response 
of  this  simple  device  structure.  The  basic  layout  of  the  de¬ 
vice  is  similar  to  the  generic  layout  illustrated  in  Fig.  13 
but  with  a  0.1  pm  gate  length.  Specihc  device  dimensions 
are  given  in  Table  2.  An  80  x  35  x  25  uniform  mesh  was 
used.  The  dc  bias  point  chosen  placed  the  active  device  well 
into  the  saturation  region  of  transistor  operation.  A  Gaussian 
pulse  excitation  with  a  3  ps  pulse  width  was  applied  over 
a  20  ps  period  to  both  gate  and  drain  electrodes  separately 
resulting  in  a  frequency  resolution  of  50  GHz.  The  source- 
gate  (input  voltage)  and  source-drain  (output  voltage)  were 


^  Springer 


J  Comput  Electron 

Table  2  Parameters  used  in  device  simulation 


drain/source  contacts  0.5  pm 

source-gate  separation  0.4  pm 

drain-gate  separation  0.4  pm 

device  length  2.0  pm 

device  width  125  pm 

active  layer  thickness  0.1  pm 

active  layer  doping  2  x  10  17  cm  3 

Schottky  barrier  height  0.8  V 

dc  gate-source  voltage  —0.5  V 

dc  drain  voltage  0.5  V 

ac  excitation  source  peak  voltage  0. 1  V 


computed  just  forward  of  the  source  plane  and  at  125  pm  re¬ 
spectively,  and  stored  at  each  time  step.  These  voltages  were 
computed  by  simply  integrating  the  v -directed  electric  field 
component  along  a  line  between  the  source-gate  and  source- 
drain  regions  as  follows: 


Fig.  19  Computed  voltage  gain  versus  frequency  for  0.1  pm  gate 
MESFET  using  a  coupled  CMC/conventional  FDTD  simulator.  A  max¬ 
imum  cutoff  frequency  of  125  GHz  is  determined  for  a  device  width  of 
125  pm 


V  =  - 


Ex  • dl , 


(17) 


where  b  —  a  represents  the  distance  between  source-gate  and 
source-drain  regions.  The  corresponding  voltage  gain  versus 
frequency  was  computed  by  taking  the  ratio  of  the  Fourier 
transforms  of  these  input  and  output  voltages  as 


Voltage  gain(&>)  =  20  log  {  z\ ^  1  (18) 

[Vin(0),Zo)\ 


where  Zi  and  z0  represent  the  distributed  and  characteristic 
impedances  and  V[n  and  V0ut  are  the  corresponding  input 
and  output  voltages.  The  result  of  this  calculation  is  given  in 
Fig.  19. 

In  addition,  the  S' -parameters  were  also  computed  over  a 
range  of  frequencies  by  taking  the  ratio  of  the  Fourier  trans¬ 
forms  of  incident  and  reflected  voltages  on  port  1  (source- 
gate)  and  port  2  (drain-source)  in  the  front  passive  section 
just  before  the  active  device  region.  Using  these  parameters, 
the  forward  current  gain,  h2\ ,  was  computed  via 


Fig.  20  Computed  current  gain  versus  frequency  for  0.1  pm  gate 
MESFET  using  a  coupled  CMC/conventional  FDTD  simulator.  An  es¬ 
timate  for  the  maximum  cutoff  frequency  of  170  GHz  is  determined 
for  a  device  width  of  125  pm 


(l-Sn)(l  +  S22)  +  Si2S2l' 

The  result  of  this  calculation  is  plotted  in  Fig.  20.  A  cutoff 
frequency  of  approximately  170  GHz  was  estimated  from 
the  curve. 

This  work  demonstrated  that  an  EM  wave  solver  could 
be  directly  coupled  to  a  full-band  particle-based  simulator 
for  the  global  modeling  of  a  microwave  transistor.  The  com¬ 
puted  voltage  gain  was  approximately  36%  larger  than  that 
predicted  by  previous  global  modeling  of  a  similar  device 


using  a  hydrodynamic-based  simulator  coupled  to  a  con¬ 
ventional  FDTD  [6,  89-92].  This  difference  is  attributed  in 
part  to  the  statistical  noise  inherent  to  particle-based  solvers 
and  also  to  the  specific  manner  in  which  the  input  and  out¬ 
put  ports  of  the  simulated  structure  were  defined.  This  latter 
issue  will  be  discussed  in  more  detail  in  the  next  section 
and  addressed  through  the  adoption  and  implementation  of 
a  typical  ground- signal-ground  (GSG)  device  layout. 

The  coupled  simulator  used  in  this  case  employed  a 
conventional  FDTD.  It  was,  therefore,  severely  limited 
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by  the  Courant  stability  requirement.  Although  the  initial 
dc  bias  point  solution  could  be  obtained  using  a  coupled 
CMC/Poisson  time  step  of  5  fs,  the  full- wave  portion  of  the 
total  simulation  was  limited  to  a  maximum  FDTD  time  step 
of  0.04  fs  or  about  two  orders  of  magnitude  smaller  than 
that  necessary  in  the  CMC  solver.  This  added  computational 
burden  has  is  resolved  with  the  implementation  of  the  ADI- 
FDTD  method,  as  described  below. 

4.3  Dual-finger  0.2  pm  gate  GaAs  MESFET 

The  2D  cross-section  of  the  3D  dual-gate  GaAs  MESFET 
simulated  using  the  CMC/ADI-FDTD  coupled  simulator  is 
shown  in  Fig.  21.  This  structure  is  simply  a  combination  of 
two  single-gate  MESFETs. 

A  141  x  25  uniform  grid  mesh  was  used  along  the  x-  and 
y -directions,  respectively.  The  choice  of  a  uniform  grid  was 
based  upon  the  fact  that  both  the  CMC  and  FDTD  solvers 
operate  on  the  same  grids  during  the  simulation.  Although 
the  CMC  grid  is  embedded  within  the  total  FDTD  grid,  its 
dimensions  and  cell  spacings  are  matched  to  those  used  in 
the  full- wave  portion  of  the  simulation.  The  use  of  a  uni¬ 
form  grid  mesh  in  the  FDTD  solver  helps  to  minimize  the 
numerical  dispersion  due  to  grid  error.  The  standard  ADI- 
FDTD  algorithm  suffers  from  the  build  up  of  both  local  and 
global  error  over  time  due  to  truncation  errors  made  in  the 
approximation  of  the  derivative  operator. 

4.3.1  Full-wave  analysis  of  the  passive  structure 

A  full-wave  simulation  of  the  passive  device  structure  was 
performed  in  order  to  generate  a  full  set  of  S -parameters  for 
comparison  with  those  obtained  during  the  coupled  simula¬ 
tions  discussed  in  the  next  section.  In  this  simulation,  the 
ADI-FDTD  algorithm  with  CPML  is  used.  The  simulation 
was  performed  using  a  141  x  60  x  35  uniform  cubic  mesh. 
These  simulations  were  run  using  an  FDTD  time  step  of  5  fs. 
A  differentiated  Gaussian  pulse  with  a  halfwidth  of  1.5  ps  is 
used  as  the  primary  excitation  source  in  the  passive  ac  simu¬ 
lation.  It  is  applied  at  3  cells  away  from  the  CPML  on  either 
the  gate  or  drain. 
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Fig.  22  Snapshots  at  (a)  2.5  ps,  (b)  5  ps,  (c)  7.5  ps,  (d)  10  ps, 
(e)  12.5  ps  and  (f)  15  ps  of  the  magnitude  of  Ex  for  the  excitation 
applied  to  the  gate 


Fig.  23  Set  of  calculated  S -parameters  for  the  passive  structure 

Snapshots  of  the  wave  propagation  resulting  from  the  ex¬ 
citation  source  being  applied  to  the  gate  are  given  in  Fig.  22. 
The  calculated  S-parameters  are  shown  in  Fig.  23. 

4.3.2  Full-band/full-wave  analysis  of  the  coupled  device 
structure 

A  fully  coupled  3D  simulation  of  the  structure  shown  in 
Fig.  21  was  performed  using  the  same  parameters  and  ex¬ 
citation  source  as  above.  At  the  start  of  the  full- wave  simu¬ 
lation,  the  CMC  device  simulation  is  allowed  to  run  coupled 
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Fig.  24  Snapshots  at  (a)  0.15  ps,  (b)  7.35  ps,  (c)  9.85  ps,  (d)  15.85  ps, 
(e)  18.85  ps  and  (f)  21.85  ps,  of  the  magnitude  of  Ex  for  the  excitation 
applied  to  the  gate 


to  the  full- wave  solver  for  5  ps  prior  to  the  application  of  the 
excitation  source  pulse.  This  allows  a  finite  amount  of  time 
for  the  resulting  ac  transients  to  propagate  throughout  the 
FDTD  grid.  These  transients  are  due  to  the  addition  of  the 
time-varying  magnetic  field  into  the  calculation  of  the  total 
Lorentz  force  in  the  CMC  portion  of  the  coupled  simulator. 
This  results  in  a  new  steady  state  (quasi-static  dc  condition). 

Snapshots  of  the  resulting  wave  propagation  are  shown  in 
Fig.  24.  At  the  beginning  of  the  simulation,  the  electric  field 
is  localized  within  the  active  region  of  the  domain  as  shown 
in  the  figures.  The  inclusion  of  the  magnetic  field  results  in  a 
redistribution  of  the  fields  out  into  the  passive  regions  of  the 
structure.  The  carrier- wave  interaction  is  also  apparent  in  the 
localized  or  granular  nature  of  the  wavefront  as  it  propagates 
through  the  active  portion  of  the  domain,  and  reflects  off  the 
open-ended  terminations  of  both  the  gate  and  drain  contacts. 

The  full  set  of  calculated  ^-parameters  are  shown  in 
Fig.  25.  The  corresponding  current  gain  /*2i,  calculated  ac¬ 
cording  to  (19),  is  plotted  in  Fig.  26.  The  unilateral  power 
gain  (UPG)  was  computed  using  the  following  expression 


UPG  = 


i  £>ii2 

|(l-|fril2)(l-|S22|2)|2 


(20) 


and  the  resulting  curve  plotted  in  Fig.  27.  The  cutoff  fre¬ 
quency  (corresponding  to  unity  current  gain)  and  /max  (cor¬ 
responding  to  unity  UPG)  are  readily  extracted  from  Figs.  26 
and  27. 


Frequency  (GHz)  Frequency  (GHz) 


Fig.  25  Set  of  calculated  5-parameters  for  the  active  structure 


Frequency  (GHz) 


Fig.  26  Current  gain  (/z2i)  versus  frequency 

In  these  simulations,  both  the  gate  and  drain  contacts 
were  assumed  to  be  simply  transmission  lines  with  charac¬ 
teristic  impedances  matched  to  that  of  the  source  and  load. 
The  CPML  used  to  truncate  the  domain  was  used  as  a  substi¬ 
tute  for  the  lumped  circuit  elements  that  would  represent  the 
source  and  load  impedances.  This  is  an  over-simplification 
of  the  problem  and  does  not  fully  account  for  the  true  im¬ 
pedances  of  the  discontinuities  present  in  actual  fabricated 
devices. 


5  Summary  and  outlook 

We  have  reviewed  recent  advances  in  the  global  model¬ 
ing  of  carrier-field  interaction  in  semiconductor  materials 
and  devices  using  coupled  EMC-FDTD  simulation  tools. 
Of  special  interest  are  applications  at  high-microwave  and 
THz  frequencies,  as  doped  semiconductor  plasma  frequen¬ 
cies  and  characteristic  carrier  scattering  rates  fall  into  this 


4?)  Springer 


J  Comput  Electron 


Fig.  27  Unilateral  power  gain  (UPG)  versus  frequency 


range.  Therefore,  common  simplifications  used  to  allow  in¬ 
dependent  consideration  of  the  particle  and  the  field  dynam¬ 
ics  are  not  justified;  rather,  efficient  modeling  tools  that  com¬ 
prehensively  couple  carrier  transport  with  full-wave  electro¬ 
dynamics  are  necessary  to  provide  an  accurate  description 
of  materials  properties. 

Coupled  EMC-FDTD  simulation  tools  combine  two 
highly  accurate  and  versatile  techniques,  and  they  present 
state-of-the-art  in  global  modeling.  The  coupled  simula¬ 
tion  was  illustrated  on  two  types  of  THz-frequency  exci¬ 
tation:  electromagnetic  irradiation  (Sect.  3)  and  ac  biasing 
(Sect.  4).  We  discussed  the  stability  and  accuracy  require¬ 
ments  for  both.  In  particular,  we  presented  simulation  data 
for  the  THz-regime  conductivity  of  doped  bulk  silicon  under 
irradiation,  and  showed  a  very  good  agreement  of  the  calcu¬ 
lated  complex  conductivity  with  experimental  data  (Sect.  3). 
In  Sect.  4.1,  ultrafast  carrier  dynamics  and  radiation  patterns 
in  thin  intrinsic  GaAs  are  faithfully  captured  via  a  CMC- 
FDTD  simulation  (CMC  being  an  outgrowth  of  ensemble 
Monte  Carlo),  while  the  same  technique  offers  a  detailed 
numerical  insight  into  the  ac  response  of  high-speed  GaAs 
MESFETs  (Sect.  4.2). 

Efficiency  of  global  modeling  tools  will  continue  to  im¬ 
prove  through  algorithm  development  and  parallelization. 
Further  developments  likely  include  analyses  of  heating  in 
high-speed  devices  by  adding  thermal  transport  to  the  cou¬ 
pled  carrier-field  simulation.  Also,  growing  interest  in  high- 
frequency  behavior  of  quantum  structures  may  stimulate 
coupling  of  quantum  transport  techniques  (Wigner-function 
simulation  or  nonequilibrium  Green’s  functions)  with  full 
wave  electrodynamics. 
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Abstract — The  experimental  results  reported  in  this  paper  sug¬ 
gest  that  single- walled  carbon  nanotubes  (SWCNTs)  have  the  po¬ 
tential  to  enhance  dielectric  contrast  between  malignant  and  nor¬ 
mal  tissue  for  microwave  detection  of  breast  cancer  and  facilitate 
selective  heating  of  malignant  tissue  for  microwave  hyperthermia 
treatment  of  breast  cancer.  In  this  study,  we  constructed  tissue- 
mimicking  materials  with  varying  concentrations  of  SWCNTs  and 
characterized  their  dielectric  properties  and  heating  response.  At 
SWCNT  concentrations  of  less  than  0.5%  by  weight,  we  observed 
significant  increases  in  the  relative  permittivity  and  effective  con¬ 
ductivity.  In  microwave  heating  experiments,  we  observed  signifi¬ 
cantly  greater  temperature  increases  in  mixtures  containing  SWC¬ 
NTs.  These  temperature  increases  scaled  linearly  with  the  effective 
conductivity  of  the  mixtures.  This  work  is  a  first  step  towards 
the  development  of  functionalized,  tumor- targeting  SWCNTs  as 
theranostic  (integrated  therapeutic  and  diagnostic)  agents  for  mi¬ 
crowave  breast  cancer  detection  and  treatment. 

Index  Terms — Breast  cancer,  carbon  nanotubes,  contrast  agent, 
dielectric  spectroscopy,  microwave  hyperthermia,  microwave 
imaging,  phantoms. 

I.  Introduction 

MICROWAVE-FREQUENCY  dielectric  contrast  between 
malignant  and  normal  tissue  in  the  breast  serves  as  the 
physical  basis  for  emerging  microwave  methods  of  detecting 
and  treating  breast  cancer.  Dielectric  contrast  leads  to  scatter- 
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ing  of  an  illuminating  microwave  signal,  which  is  exploited  for 
breast  cancer  detection,  and  selective  absorption  of  incident  mi¬ 
crowave  power,  which  is  exploited  for  localized  hyperthermia 
treatment.  The  effective  dielectric  properties  of  breast  tissue 
are  influenced  at  microwave  frequencies  by  endogenous  polar 
molecules,  such  as  free  and  bound  water,  peptides,  and  proteins. 
Consequently,  the  dielectric  properties  depend  on  the  type  and 
physiological  state  of  the  tissue.  A  recent  study  conducted  by  the 
University  of  Wisconsin-Madison  and  the  University  of  Calgary 
showed  that  the  endogenous  dielectric  contrast  between  malig¬ 
nant  and  normal  adipose-dominated  tissues  in  the  breast  is  as 
large  as  10:1,  while  the  contrast  between  malignant  and  normal 
glandular/fibroconnective  tissues  is  no  more  than  about  10%  [1]. 

The  effective  dielectric  properties — both  the  dielectric  con¬ 
stant  and  effective  conductivity — can  also  be  influenced  by  ex¬ 
ogenous  molecules  introduced  as  contrast  agents.  For  example, 
another  recent  study  found  that  tissue-mimicking  (TM)  mix¬ 
tures  containing  encapsulated  microbubbles  at  clinically  rele¬ 
vant  concentrations,  similar  to  those  used  in  ultrasound  imag¬ 
ing,  exhibit  significantly  lower  dielectric  properties  than  pure 
TM  materials  [2].  Here,  we  propose  the  use  of  single- walled 
carbon  nanotubes  (SWCNTs)  as  a  diagnostic  and  therapeutic 
agent — an  integrated  “theranostic”  nanoplatform  [3] — for  both 
microwave  detection  and  treatment  of  breast  cancer.  The  size 
and  unique  physiochemical  properties  of  SWCNTs  make  them 
ideal  candidates  for  tumor-targeting  applications  [4].  SWCNT- 
based  contrast  agents  have  already  shown  promise  for  molecular 
imaging  by  magnetic  resonance,  positron  emission  tomography, 
nuclear,  and  photoacoustic  imaging  modalities  [5]-[7]. 

Our  present  investigation  is  motivated  by  the  hypothesis  that 
the  accumulation  of  biocompatible  carbon  nanotubes  in  tumors 
will  significantly  increase  the  microwave-frequency  dielectric 
properties  of  the  tumor.  Exogenous  particles  passively  accumu¬ 
late  via  an  enhanced  permeability  and  retention  effect  in  tumor 
vasculature  [8].  Functionalization  of  these  particles  with  tumor¬ 
targeting  biomarkers  may  amplify  their  accumulation  properties 
in  tumors  [9].  The  presence  of  carbon  nanotubes  in  or  near  a 
malignant  lesion  will  result  in  changes  to  the  dielectric  prop¬ 
erties  of  malignant  tissue.  If  these  changes  are  significant,  they 
can  be  used  to  enhance  the  sensitivity  of  low-power  microwave 
imaging  of  breast  cancer  via  differential  imaging  [10],  [11]  and 
to  enhance  the  selectivity  of  higher  power  microwave  thermal 
therapy. 
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In  this  paper,  we  report  the  results  of  wideband  (0.6-20  GHz) 
dielectric  spectroscopy  measurements  and  heating  efficiency 
experiments  (at  3  GHz)  using  TM  materials  mixed  with  var¬ 
ious  concentrations  of  SWCNTs  that  have  been  identified  to 
be  non-toxic  in  mice  [12].  Our  study  addresses  the  dielectric 
and  thermal  effect  of  SWCNTs  over  the  frequency  range  that 
is  of  interest  for  microwave  imaging  and  hyperthermia  treat¬ 
ment  of  breast  cancer.  Frequencies  between  0.5  and  3  GHz  are 
commonly  employed  in  microwave  imaging  via  inverse  scatter¬ 
ing  [13]  while  the  ultra- wideband  range  of  3.1-10.6  GHz  is  of 
interest  in  tissue-penetrating  radar  imaging  [14].  Microwave- 
induced  thermoacoustic  imaging  has  been  explored  using 
434  MHz  [15]  and  3  GHz  [16].  Microwave  hyperthermia  treat¬ 
ment  of  breast  cancer  has  been  conducted  at  915  MHz  [17], 
although  slightly  higher  frequencies  (1. 5-4.0  GHz)  have  been 
shown  recently  to  be  optimal  for  tightly  focusing  microwave 
energy  in  the  breast  [18].  We  choose  3  GHz  for  our  heating  ex¬ 
periments  because  of  its  relevance  to  both  microwave-induced 
thermoacoustic  imaging  and  hyperthermia  treatment.  Our  work 
complements  previous  radiofrequency  heating  studies  [19]  as 
well  as  previous  microwave-frequency  characterizations  of  the 
dielectric  properties  of  mixtures  of  SWCNTs  or  multi  walled 
carbon  nanotubes  in  various  host  media  of  interest  for  electro¬ 
magnetic  shielding  and  radar  absorption  [20]-[26]. 


II.  Materials  and  Methods 

A.  Construction  of  Samples 

The  TM  materials  were  constructed  from  oil-in-gelatin  dis¬ 
persions  following  the  procedure  outlined  in  [27] .  The  dielectric 
properties  of  these  materials  can  be  customized  to  mimic  the 
properties  of  a  variety  of  human  soft  tissues  by  controlling  the 
concentrations  of  gelatin,  safflower  oil,  kerosene,  and  preserva¬ 
tives.  The  results  reported  in  [1]  show  that  a  10%-oil  mixture 
adequately  replicates  the  microwave  properties  of  malignant 
breast  tissue  over  our  frequency  range  of  interest  (0.6  GHz  to 
20  GHz).  Furthermore,  these  materials  are  relatively  inexpen¬ 
sive  to  fabricate  and  possess  long  term  stability. 

For  this  study  we  constructed  10%-oil  TM  materials  with 
various  concentrations  of  SWCNTs.  The  SWCNTs  were 
synthesized  by  a  commercial  vendor  (Cheap  Tubes  Inc., 
Brattleboro,  VT)  using  a  chemical  vapor  deposition  technique 
and  acid  treated  for  purification.  According  to  the  manufac¬ 
turer’s  specifications,  the  nanotubes  were  1-2  nm  in  diameter 
and  5-30  pm  in  length,  and  were  composed  of  mainly  SWCNTs 
(>90%  by  weight)  with  trace  amounts  of  multi  walled  carbon 
nanotubes  and  metal  catalyst  (e.g.,  cobalt,  1%  by  weight).  We 
mixed  various  concentrations  (1,2,  and  3  mg/mL)  of  SWCNTs 
into  a  1%  by  weight  mixture  of  Pluronic  (FI 27)  and  deionized 
water  using  a  probe  sonicator  for  25  min  at  a  120  W  power 
level.  Then,  we  substituted  these  mixtures  for  the  pure  deion¬ 
ized  water  in  the  TM  recipe  described  in  [27].  The  resulting 
percent-by- weight  concentrations  of  SWCNTs  in  the  samples 
were  0.07%,  0.15%,  and  0.22%,  respectively.  For  reference,  we 
also  constructed  a  pure  TM  mixture  with  the  same  amount  of 
Pluronic  as  the  other  samples,  but  with  no  SWCNTs. 


Tempera  Lure 


Fig.  1.  Experimental  setup  for  characterizing  the  heating  response  of  a  TM 
sample. 

B.  Characterization  of  Dielectric  Properties 

The  dielectric  properties  of  the  TM  samples  were  character¬ 
ized  using  a  well  established  open-ended  coaxial  probe  tech¬ 
nique  described  in  [28].  The  measurement  uncertainty  of  this 
technique  is  no  more  than  approximately  10%  [28].  We  poured 
the  liquid  oil-gelatin  mixtures,  described  in  Section  II-A,  into 
small  covered  glass  jars  and  allowed  the  mixture  to  gel  un¬ 
perturbed.  We  positioned  the  tip  of  the  precision  coaxial  probe 
against  the  surface  of  the  sample  [27] .  We  conducted  these  mea¬ 
surements  at  room  temperature  (approximately  22  °C)  on  three 
TM  samples  of  each  SWCNT  concentration. 

C.  Characterization  of  Heating  Response 

The  configuration  for  the  heating  experiment  is  shown  in 
Fig.  1.  We  prepared  the  samples  for  the  heating  experiments  by 
pouring  the  liquid  TM  mixture,  as  described  in  Section  II-A, 
into  a  1.1-mm-inner-diameter  glass  capillary  tube.  The  mixture 
was  allowed  to  gel  around  a  fiber-optic  temperature  probe  con¬ 
nected  to  a  calibrated  fluoroptic  thermometer  (Luxtron,  3100). 
The  capillary  tube  was  inserted  through  a  small  hole  drilled  into 
the  center  of  the  broad  wall  of  an  S-band  (WR-284)  rectangular 
waveguide  with  cross-sectional  dimensions  of  72  mm  x  34  mm. 
The  vertical  extent  of  the  sample  in  capillary  tube  spanned 
the  entire  height  of  the  waveguide.  A  microwave  synthesizer 
(Agilent,  83623B)  and  amplifier  (Mini-Circuits,  ZHL-42W) 
generated  1  W  of  continuous  microwave  power  at  3  GHz  that 
was  delivered  to  the  sample  via  the  waveguide.  The  waveguide 
was  terminated  with  a  matched  load  to  suppress  standing  waves 
and  ensure  single-pass  electromagnetic  wave  illumination.  Net¬ 
work  analyzer  measurements  verified  that  the  small  hole  drilled 
in  the  waveguide  and  the  presence  of  the  sample  did  not  sig¬ 
nificantly  perturb  the  fields  inside  the  waveguide;  namely,  these 
waveguide  alterations  increased  |5n  |  by  approximately  0.6  dB. 
The  fluoroptic  thermometer  recorded  the  time- varying  temper¬ 
ature  of  each  sample  as  it  was  heated.  We  heated  each  sample 
for  3  min,  then  turned  off  the  microwave  source  and  allowed  it 
to  cool  for  5  min. 

III.  Results  and  Discussion 
A.  Dielectric  Properties 

Fig.  2(a)  and  (b)  shows  the  relative  permittivity  and  effective 
conductivity  of  the  TM  materials  with  various  concentrations  of 
SWCNTs  over  a  frequency  range  of  0.6-20  GHz.  Each  wide¬ 
band  curve  represents  the  average  properties  (relative  permittiv¬ 
ity  or  effective  conductivity)  of  the  three  samples  with  the  same 
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Fig.  2.  (a)  Relative  permittivity  and  (b)  effective  conductivity  of  TM  materials 
with  varying  concentrations  of  SWCNTs,  measured  from  0.6  to  20  GHz.  Each 
curve  represents  the  average  properties  of  three  different  samples  of  the  same 
concentration.  The  bars  span  the  maximum  and  minimum  values  at  specific 
frequencies.  Both  the  permittivity  and  effective  conductivity  of  the  TM  materials 
increase  with  increasing  concentration  of  SWCNTs. 


TABLE  I 

Summary  of  Data  Presented  in  Fig.  2  at  3  GHz 


Percent  by 
weight 

concentration 
of  SWCNTs 

Average 
er  (3  GHz) 

Average 
a  (3  GHz) 

Average 
percent 
change 
er  (3  GHz) 

Average 

percent 

change 

CT  (3  GHz) 

None 

53.3 

2.3  S/m 

- 

- 

0.07  % 

57.9 

2.7  S/m 

9  % 

16  % 

0.15  % 

63.7 

3.2  S/m 

19  % 

38  % 

0.22  % 

73.0 

4.2  S/m 

37  % 

81  % 

SWCNT  concentration.  The  vertical  bars  span  the  maximum 
and  minimum  values  at  specific  frequencies.  Table  I  summa¬ 
rizes  the  dielectric  properties  of  the  different  TM  mixtures  at 
3  GHz.  This  particular  frequency  is  of  interest  for  microwave 
inverse  scattering  [10],  [11],  microwave-induced  thermoacous¬ 
tic  imaging  [16],  and  microwave  hyperthermia  [18].  Both  Fig.  2 
and  Table  I  show  that  the  permittivity  and  conductivity  increase 
with  small  increases  in  the  concentration  of  SWCNTs. 

These  substantial  changes  in  the  dielectric  properties  caused 
by  the  carbon  nanotubes  will  have  a  significant  impact  on  mi¬ 
crowave  scattering  and  absorption.  The  results  reported  here 
were  used  in  a  microwave  imaging  computational  test  bed  in  or¬ 
der  to  determine  the  impact  of  these  changes  in  microwave  breast 
cancer  detection  [10],  [11].  Briefly,  microwave  inverse  scatter¬ 
ing  techniques  were  used  to  reconstruct  images  of  anatomically 
realistic  numerical  breast  phantoms  with  a  compact  malignant 


Fig.  3.  Microwave  heating  response  of  TM  materials  with  various  concentra¬ 
tions  of  SWCNTs.  Each  curve  shows  the  temperature  profile  of  a  sample  that 
was  heated  via  3-GHz  microwave  illumination  for  3  min  and  allowed  to  cool  for 
5  min.  The  maximum  temperature  of  the  TM  mixtures  increases  with  increasing 
concentration  of  SWCNTs. 


TABLE  II 

Summary  of  Data  Presented  in  Fig.  3 


Percent  by 
weight 

concentration 
of  SWCNTs 

Average 

steady-state 

temperature 

Average 

temperature 

rise 

AT  increase 
compared  to 
SWCNT-free 
sample 

None 

30.9°C 

6.9°C 

- 

0.07  % 

33.0°C 

9,0°C 

2.1°C 

0.15  % 

34.5°C 

10.5°C 

3.6°C 

0.22  % 

36.9°C 

12.9°C 

6.0°C 

lesion.  Two  images  were  reconstructed:  an  image  of  a  phantom 
with  the  endogenous  dielectric  properties  of  the  tumor  and  an¬ 
other  with  elevated  dielectric  properties  values  due  to  the  carbon 
nanotubes.  A  differential  image  was  formed  by  subtracting  the 
two  images  obtained  from  pre-  and  postcontrast  measurements. 
Using  this  differential  imaging  technique,  Shea  et  al.  [10],  [11] 
reported  detection  of  previously  undetectable  tumors  in  four  dif¬ 
ferent  classes  of  numerical  phantoms  that  ranged  from  “mostly 
fatty”  to  “extremely  dense”  in  radiographic  density.  These  re¬ 
sults  suggest  that  the  changes  in  dielectric  properties  shown  in 
Fig.  2  are  significant  enough  to  dramatically  improve  the  sensi¬ 
tivity  of  microwave  imaging. 

B.  Heating  Response 

Fig.  3  shows  the  time-dependent  temperature  rise  of  the  var¬ 
ious  TM  mixtures  under  continuous  microwave  illumination. 
The  individual  temperature  value  on  each  curve  was  a  rolling 
average  of  five  different  temperature  measurements  of  a  single 
sample  that  was  recorded  every  0.25  s.  For  each  TM  mixture 
concentration,  measurements  were  performed  on  three  different 
samples  of  each  TM  mixture.  Fig.  3  demonstrates  the  mea¬ 
surement  repeatability  for  each  concentration,  as  well  as  the 
statistical  significance  of  the  differences  observed  across  con¬ 
centrations.  Table  II  summarizes  the  principal  characteristics  of 
the  results  shown  in  Fig.  3.  These  data  show  that  as  the  concen¬ 
tration  of  SWCNTs  in  the  TM  material  increases,  the  maximum 
temperature  of  that  sample  also  increases.  This  increase  in  tem¬ 
perature  for  each  TM  mixture  scales  roughly  linearly  with  the 
increase  in  conductivity  reported  in  Table  I,  as  expected  for  a 
configuration  in  which  various  samples  are  exposed  to  approxi¬ 
mately  the  same  electric  field  intensity.  In  other  experiments,  we 
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have  observed  that  microwave-heated  capillary  tubes  of  liquids 
and  SWCNT-containing  liquids  similarly  exhibit  a  temperature 
rise  proportional  to  their  effective  conductivity.  Therefore,  this 
heating  response  scaling  appears  to  be  a  general  trend,  not  re¬ 
stricted  to  just  TM  material  mixtures. 

The  enhanced  heating  shown  in  Fig.  3  due  to  the  carbon 
nanotubes  is  of  interest  to  microwave-induced  thermoacoustic 
imaging  and  microwave  hyperthermia  of  breast  cancer.  Our  re¬ 
sults  suggest  that  the  increase  in  heating  response  will  lead  to 
a  larger  thermoacoustic  response  for  a  malignant  tumor  that 
is  infused  with  carbon  nanotubes  in  comparison  to  a  tumor  in 
a  precontrast  stage.  In  addition,  the  results  of  Fig.  3  suggest 
that  malignant  tissue  that  has  preferentially  accumulated  carbon 
nanotubes  will  heat  more  efficiently;  this  selective  heating  is  a 
desired  property  for  microwave  hyperthermia  of  breast  cancer. 
As  seen  in  Fig.  2(b),  various  samples  exhibit  a  large  spread  in 
effective  conductivity  over  the  range  of  0.6-20  GHz.  Therefore, 
we  expect  a  similar  heating  response  to  those  seen  in  Fig.  3  for 
these  samples  over  the  0.6-20  GHz  frequency  range. 

IV.  Conclusion 

We  characterized  the  dielectric  properties  of  TM  materials 
with  several  different  concentrations  of  SWCNTs.  Our  results 
indicate  that  low  concentrations  of  SWCNTs  significantly  im¬ 
pact  the  dielectric  properties  and  heating  response  of  the  TM 
materials.  For  example,  at  3  GHz,  SWCNTs  concentrations  as 
small  as  0.22%  by  weight  increased  the  relative  permittivity  of 
the  TM  material  by  37%  and  the  effective  conductivity  by  81%. 
In  our  heating  experiments,  this  concentration  of  SWCNTs  led 
to  an  average  steady-state  temperature  rise  that  was  6  °C  higher 
than  the  rise  observed  in  the  TM  material  without  SWCNTs. 
These  results  suggest  that  SWCNTs  may  enhance  contrast  for 
microwave  imaging  and  facilitate  selective  microwave  heating 
for  treatment  of  breast  cancer. 
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Abstract 

The  physical  basis  for  microwave  breast  cancer  detection  is  the  dielectric-properties  contrast  between 
malignant  and  normal  breast  tissue.  The  Wisconsin-Calgary  study  showed  that  this  contrast  is  as  high  as  10:1  in 
fatty  breast  tissue  but  no  more  than  10%  in  fibroglandular  tissue.  We  are  investigating  the  feasibility  of  air-filled 
microbubbles  as  contrast  agents  for  enhancing  the  malignant-to-normal  dielectric  contrast.  Our  initial  numerical 
studies  suggest  that  the  presence  of  moderate  volume  fractions  of  microbubbles  reduces  the  effective  dielectric 
properties  of  the  tumor  by  as  much  as  30%,  potentially  improving  detection  efficacy  via  differential  imaging. 

1.  Introduction 

Emerging  microwave  breast  cancer  detection  techniques  seek  to  exploit  a  large  dielectric-properties 
contrast  between  malignant  and  normal  breast  tissues.  Our  recently- completed  large-scale  dielectric 
characterization  study  of  human  breast  tissue  samples  conducted  over  the  0.5-20  GHz  frequency  range 
demonstrated  that  the  dielectric-properties  contrast  between  malignant  and  normal  adipose-dominated  breast 
tissues  is  as  large  as  10:1,  while  the  dielectric-properties  contrast  between  malignant  and  normal  fibroglandular 
tissues  is  no  more  than  approximately  10%  [1].  Since  the  vast  majority  of  breast  tumors  originate  within 
fibroglandular  breast  tissue,  our  results  suggest  that  the  malignant  lesion  is  a  weakly  scattering  target  within  a 
high  clutter  environment.  To  address  this  challenging  imaging  scenario,  we  are  assessing  the  feasibility  of  using 
dielectric  or  conducting  micro-  and  nano-particles  as  contrast  agents  for  enhancing  the  dielectric-properties 
contrast  between  the  tumor  and  surrounding  normal  fibroglandular  tissue.  The  specific  particles  we  have  chosen 
for  this  analysis  are  air-filled  microbubbles,  such  as  the  ones  used  as  ultrasound  contrast  agents,  and  metallic 
nanoparticles.  Once  introduced  into  the  breast,  the  particles  are  assumed  to  accumulate  in  the  malignancy  via 
passive  and/or  active  mechanisms  [2].  The  efficacy  of  tumor  detection  is  expected  to  be  enhanced  via 
differential  imaging  of  the  breast.  In  this  paper  we  present  the  results  of  our  initial  numerical  studies  of 
dielectric-properties  contrast  enhancement  for  the  case  of  air- filled  microbubble  contrast  agent  particles.  A 
comparison  between  experimental  and  numerical  results  is  also  presented. 

2.  Methodology 

The  accumulation  of  contrast  agent  particles  in  a  breast  tumor  is  modeled  as  a  binary  mixture  consisting 
of  a  certain  volume  fraction  (VF)  of  inclusions  (microbubbles)  in  a  background  medium  (malignant  breast 
tissue).  We  implemented  both  a  2D  quasi-static  finite- element  method  (QS-FEM),  as  well  as  a  2D  finite- 
difference  time-domain  (FDTD)  approach.  The  advantages  of  the  former  method  are  that  it  is  simple  and  fast. 
However,  it  does  not  account  for  the  physical  size  of  the  inclusions,  but  simply  assumes  that  the  particles  are 
much  smaller  than  the  wavelength  of  the  electromagnetic  wave  in  the  mixture.  Since  size-related  effects,  such  as 
skin  depth,  are  particularly  important  for  finite-conductivity  inclusions,  such  as  metallic  nanoparticles,  we 
conducted  full- wave  FDTD  simulations  to  validate  the  QS-FEM  technique  and  obtain  supplemental  information. 
Here  we  address  the  simpler  case  of  air-filled  microbubble  inclusions,  since  the  results  for  this  case  are  not 
confounded  by  skin  depth  and  percolation  effects. 

For  both  numerical  approaches,  the  dielectric  properties  of  the  background  medium  in  the  computational 
domain  [Figures  1(a)  and  2(a)]  were  set  to  those  of  malignant  breast  tissue  at  5  GHz  (£r^  =50.65,  (J^=4.84 

S/m)  [1].  The  dielectric  properties  of  the  air- filled  microbubbles  were  set  to  those  of  free  space  (£ri=l,  Gi  =0 

S/m).  The  randomly-positioned  inclusions  were  allowed  to  touch,  but  not  overlap.  For  simplicity,  the  particles 
did  not  intersect  the  edges  of  the  computational  domain.  For  every  YF  of  inclusions  (5,  10,  and  20%),  100 
geometries  with  random  particle  locations  were  generated  and  simulated  in  order  to  capture  the  effects  of 
different  inclusion  distributions  in  the  ensemble  average  of  effective  dielectric  properties  values.  For  every  VF, 
an  average  effective  dielectric  constant  and  conductivity  was  computed. 
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2.1.  QS-FEM  approach 

The  2D  QS-FEM  approach,  based  on  [3],  was  implemented  using  the  partial  differential  equation 
toolbox  in  Matlab.  The  QS-FEM  solved  Laplace’s  equation  for  the  geometry  and  boundary  conditions  of  Figure 
1(a).  Once  the  electrostatic  potential  0(r)  was  computed  at  every  point  r  in  the  computational  domain,  the 
stored  and  dissipated  energies  in  the  capacitor  of  Figure  1(a)  (henceforth  referred  to  as  a  complex  stored  energy) 
were  found  using  W*  =  1/2£q  j*£* (r)[V0(r)]2d2r  ,  where £* (r)  =  ^(r)  —  j <j(y)I(cQ£q)  is  the  complex 
permittivity  at  every  point  in  the  computational  domain,  £*0  is  the  permittivity  of  free  space,  CO  is  the  angular 
frequency,  and  j  =  4— 1  .  Assuming  that  the  complex  stored  energy  of  the  capacitor  in  Figure  1(a)  is  equivalent 
to  that  of  Figure  1(b),  the  effective  dielectric  properties  in  1(b)  were  found  using  £*eff  =  2 W*j ) . 


The  top  plate  of  the  capacitor  in  Figure  1(a)  was  held  at  a  constant  potential  of  0O=1  V,  while  the 
bottom  plate  was  held  at  0  V.  In  order  to  avoid  fringing  fields  at  the  edges  of  the  structure,  Neumann  boundary 
conditions  were  enforced  on  the  left  and  right  sides  of  the  computational  domain  [Figure  1(a)].  The  size  of  the 
computational  domain  was  100  by  100  units,  while  the  radius  of  the  inclusions  was  2  units  (the  actual  physical 
dimensions  are  irrelevant).  A  mesh  conformal  to  the  geometry  was  automatically  generated  by  the  solver. 


dn 


(a) 


5*=Eefl 

^5T 


(b) 


Figure  1.  (a)  Diagram  of  the  geometry  simulated  using  QS-FEM.  (b)  Effective  medium  formulation.  8b: 
background  dielectric  properties,  £*:  inclusion  dielectric  properties,  £eff:  effective  dielectric  properties. 


2.2.  FDTD  approach 

The  2D  TEZ  FDTD  formulation  used  to  simulate  the  TEM  parallel-plate  waveguide  of  Figure  2(a)  was 
similar  to  that  of  [4].  The  computed  field  components  were  Ex,  Ey,  and  Hz.  This  polarization  was  chosen  to 
match  the  polarization  of  the  QS-FEM  simulations.  The  FDTD-computed  Ey  was  recorded  at  10  equally- spaced 
grid  cells  along  the  y-dimension  of  the  waveguide  at  observation  planes  1  and  2  [see  Figure  2(a)].  Next,  we 
calculated  the  attenuation  ( a s  )  and  propagation  ( jBs  )  constants  (assuming  that  the  TEM  mode  in  the  mixture 
“slab”  was  not  perturbed)  using  the  FDTD-computed  Fourier  transforms  of  Ey  at  planes  1  and  2.  Finally,  we 
extracted  the  effective  dielectric  properties  of  the  mixture  from  the  attenuation  and  propagation  constants.  For 
every  simulation,  the  effective  dielectric  properties  calculated  from  data  at  the  ten  observation  points  along  the 
y-dimension  were  averaged. 

We  used  a  modulated  Gaussian  source,  with  a  center  frequency  of  5  GHz,  and  a  1/e  width  of 
approximately  318  ns.  Every  simulation  was  run  as  long  as  necessary  to  sufficiently  dissipate  the  energy  inside 
the  grid.  The  TEM  mode  was  excited  within  the  waveguide  by  assigning  a  current  source  to  the  Ey  components 
along  the  source  plane  [see  Figure  2(a)]. 

The  plate  separation  was  Ny=50  grid  cells.  The  mixture  “slab”  of  width  Nx=100  grid  cells  was  placed  in 
the  center  of  the  waveguide,  where  L1=L2=50  grid  cells.  The  waveguide  was  terminated  with  perfectly  matched 
layer  absorbing  boundary  conditions.  The  top  and  bottom  walls  of  the  waveguide  were  assumed  to  be  perfect 
electric  conductors.  Since  the  microbubbles  used  as  ultrasound  contrast  agents  are  approximately  1-5  pm  in 
diameter,  we  set  the  diameter  of  the  particles  in  the  FDTD  simulations  to  be  5  pm.  In  order  to  properly  resolve 
the  particles,  the  grid  cell  size  was  1.0  pm.  We  used  the  Yu-Mittra  technique  [5]  to  compute  the  local  dielectric 
properties  for  the  electric-field  components  that  cross  the  edges  of  the  inclusions. 
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Figure  2.  (a)  Diagram  of  the  parallel-plate  waveguide  simulated  using  FDTD.  (b)  Effective  medium  formulation. 

2.3.  Comparison  with  experimental  studies 

In  addition  to  comparing  the  QS-FEM  numerical  results  with  FDTD,  we  compared  our  numerical  results 
with  experimental  data  for  the  cases  of  5%  and  10%  by  weight  mixtures  of  air- filled  microbubbles  in  ethylene 
glycol  [6].  For  this  comparison,  the  dielectric  properties  of  the  background  medium  in  the  QS-FEM  simulations 
were  set  to  those  of  ethylene  glycol  at  5  GHz.  The  density  of  the  glass-shell,  air-filled  microbubbles  (average 
diameter:  18  jam)  is  0.6  g/mL,  so  the  VF  of  microbubbles  in  the  5%  and  10%  by  weight  mixtures  were 
approximately  8.3%  and  16.7%,  respectively.  Ethylene  glycol  was  chosen  as  an  inexpensive,  readily  available 
reference  liquid,  with  dielectric  properties  representative  of  those  of  biological  tissues.  Immediately  after 
mixing,  the  dielectric  properties  of  the  mixture  were  measured  from  0.5  to  20  GHz  using  an  open-ended  coaxial 
probe  [6].  The  measurements  were  repeated  five  times  for  three  different  probe  positions  within  the 
measurement  beakers,  and  the  dielectric  properties  recorded  for  all  trials  were  averaged. 

3.  Comparison  between  QS-FEM  and  FDTD  simulations 

Table  1  shows  the  numerically  computed  effective  dielectric  properties,  averaged  over  100  iterations,  of 
three  different  mixtures  (VF=5%,  10%,  and  20%)  of  microbubbles  in  a  background  medium  representing 
malignant  breast  tissue.  For  every  VF  of  inclusions,  the  variability  in  the  effective  dielectric  properties  across 
iterations  is  very  small,  no  more  than  about  7%.  In  addition,  we  found  excellent  agreement  between  the  QS- 
FEM-computed  and  FDTD-computed  effective  dielectric  properties  of  the  mixtures.  In  fact,  for  all  three  VFs 
investigated  in  this  study,  the  QS-FEM-computed  effective  dielectric  properties  were  within  approximately  1% 
of  the  FDTD-computed  effective  dielectric  properties,  indicating  that  the  simpler  QS-FEM  formulation  is  very 
accurate  for  modeling  the  dielectric-properties  changes  due  to  microbubble  inclusions.  Furthermore,  we  found 
that  moderate  VFs  of  microbubbles  have  noticeable  effects  on  the  effective  dielectric  properties  of  malignant 
breast  tissue  (Table  1).  For  example,  a  20%  VF  of  microbubbles  resulted  in  an  approximately  30%  reduction  in 
the  effective  dielectric  properties  of  the  tissue. 


Table  1.  Comparison  between  QS-FEM-computed  and  FDTD-computed  effective  dielectric  properties  of 
mixtures  of  air-filled  microbubbles  in  breast  tumor. 


Effective  dielectric  constant  at  5  GHz 

Effective  conductivity  at  5  GHz  (S/m) 

VF  (%) 

QS-FEM 

FDTD 

%  change 

QS-FEM 

FDTD 

%  change 

0 

50.65 

50.65 

- 

4.84 

4.84 

- 

5 

46.09 

45.87 

-9 

4.39 

4.36 

-10 

10 

41.76 

41.81 

-17 

3.96 

3.96 

-18 

20 

34.46 

34.16 

-32 

3.23 

3.22 

-33 

4.  Comparison  between  numerical  and  experimental  studies 

Figure  3  shows  the  effective  dielectric  properties  of  a  5%  and  10%  by  weight  mixture  of  microbubbles  in 
ethylene  glycol.  The  solid  and  dotted  lines  represent  experimental  data,  while  the  symbols  represent  QS-FEM- 
computed  effective  dielectric  properties  at  5,  10,  and  15  GHz.  We  found  excellent  agreement  between  the 
experimental  and  numerical  results  -  the  QS-FEM-computed  results  are  within  about  8.5%  of  the  experimental 
results.  This  small  discrepancy  can  be  attributed  to  the  simplified,  2D  nature  of  the  QS-FEM  simulations  and  the 
uncertainties  in  the  actual  VF  of  inclusions  in  the  mixtures  and  the  exact  dielectric  properties  of  the 
microbubbles.  The  fact  that  the  actual  dielectric  properties  of  the  microbubbles  are  slightly  higher  than  the 
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assumed  value  of  free  space  and  are  spherical  rather  than  cylindrical  may  explain  why  the  QS-FEM-predicted 
decrease  in  the  effective  dielectric  properties  is  somewhat  larger  than  the  experimentally  measured  decrease. 


Frequency  (GHz)  Frequency  (GHz) 

Figure  3.  Experimentally- measured  and  QS-FEM  computed  effective  (a)  dielectric  constant  and  (b)  conductivity 
of  mixtures  consisting  of  5%  and  10%  by  weight  air- filled  microbubbles  in  ethylene  glycol. 

5.  Conclusions 

We  have  conducted  numerical  studies  to  investigate  the  effects  of  air-filled  microbubbles  on  the 
dielectric  properties  of  malignant  breast  tissue.  We  demonstrated  that  the  QS-FEM  technique  accurately  models 
the  change  in  dielectric  properties  of  a  dielectric  medium  representing  malignant  breast  tissue  due  to  the 
presence  of  microbubbles.  We  found  that  a  5-20%  VF  of  microbubbles  can  lower  the  dielectric  properties  of 
malignant  tissue  by  as  much  as  30%.  Through  another  research  effort  in  our  group,  we  have  explored  the 
impact  of  this  dielectric  properties  reduction  on  microwave  images  generated  for  3D  anatomically  realistic 
numerical  breast  phantoms  using  an  inverse  scattering  algorithm.  Results  from  those  studies  demonstrate  the 
feasibility  of  tumor  detection  via  differential  (pre-  and  post-contrast-agent)  imaging. 
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Microwave-frequency  dielectric  contrast  between  malignant  and  normal  tissue  in  the  breast 
serves  as  the  physical  basis  for  emerging  microwave  methods  of  detecting  and  treating  breast 
cancer.  Dielectric  contrast  leads  to  scattering  of  an  illuminating  microwave  signal,  which  is 
exploited  for  breast  cancer  detection,  and  selective  absorption  of  incident  microwave  power, 
which  is  exploited  for  localized  hyperthermia  treatment.  The  endogenous  dielectric  properties 
of  breast  tissue  are  dependent  on  the  physiological  type  and  state  of  the  tissue  (Lazebnik,  et 
al.,  Phys.  Med.  Bio.,  vol.  52,  pp.  2637-2656  and  6093-6115,  2007).  The  effective  dielectric 
properties  may  also  be  influenced  by  exogenous  molecules  introduced  as  contrast  agents.  We 
propose  the  use  of  single-walled  carbon  nanotubes  (SWCNTs)  as  a  diagnostic  and  therapeutic 
agent  -  an  integrated  “theranostic”  nanoplatform  -  for  both  microwave  detection  and 
treatment  of  breast  cancer  (Mashal,  et  al.,  IEEE  TBME  Letters ,  accepted). 

Here,  we  report  the  results  of  wideband  (0.5-10  GHz)  dielectric  spectroscopy  measurements 
and  heating  efficiency  experiments  (0.5-1  GHz)  of  SWCNT-water  dispersions.  Our  study 
addresses  the  dielectric  and  thermal  effect  of  biomedically  relevant  concentrations  of 
SWCNTs  over  the  frequency  range  that  is  of  interest  for  microwave  imaging  and  hyperthermia 
treatment  of  breast  cancer. 

We  created  the  SWCNT  dispersions  by  sonicating  SWCNTs  in  a  mixture  of  water  and 
polyethylene  glycol  (PEG).  FDA-approved  PEG  is  used  to  non-covalently  functionalize  the 
SWCNTs,  thereby  ensuring  that  they  remain  biocompatible  and  well  dispersed  with  no 
aggregation.  We  measured  the  dielectric  properties  using  a  precision  open-ended  coaxial 
probe  technique  (Popovic,  et  al.,  IEEE  Trans.  Microwave  Theory  Tech.,  vol.  53,  pp.  1713- 
1722,  2005).  For  a  concentration  of  2  mg/mL  of  SWCNTs,  we  saw  an  increase  in  permittivity 
of  about  10%  and  an  increase  in  conductivity  of  about  25%  when  compared  to  a  pure  DSPE- 
PEG  solution  across  the  0.5  GHz  to  3  GHz  frequency  range.  Next  we  measured  the  impact  of 
the  SWCNTs  on  the  heating  response  of  a  water-based  mixture.  A  sample  of  the  SWCNT 
dispersion  was  introduced  into  a  1 . 1 -mm-inner-diameter  glass  capillary  tube.  The  capillary 
tube  was  inserted  through  a  small  hole  drilled  into  the  center  of  a  hardline  coax  and  was 
illuminated  with  continuous  microwave  energy  at  discrete  frequencies  between  0.5  GHz  to  1 
GHz.  A  fluoroptic  thermometer  was  used  to  record  the  temperature  during  each  experiment. 
We  observed  the  maximum  final  temperature  increase  by  a  factor  of  three  at  0.5  GHz  and  a 
factor  of  two  at  1  GHz.  This  increase  in  temperature  scales  approximately  with  the  changes 
observed  in  the  effective  conductivities  at  these  frequencies. 

The  enhanced  dielectric  properties  and  heating  response  suggest  that  SWCNTs  are  promising 
theranostic  agents  for  microwave  detection  and  thermal  treatment  of  breast  cancer.  We  have 
previously  shown  that  the  CNT-enhanced  dielectric  contrast  greatly  improves  the  sensitivity  of 
microwave  breast  imaging  (Shea,  et  al.,  ANTEM/URSI  2009).  Here,  we  will  highlight  their 
impact  on  hyperthermia  performance  using  computational  electromagnetic  and  thermal 
simulations  of  numerical  breast  phantoms  containing  tumors  targeted  with  CNTs. 
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Semiconductor  material  properties  at  terahertz 
(THz)  and  sub-THz  frequencies  cannot  be  well  de¬ 
scribed  by  the  models  used  at  low  frequencies.  Bulk 
doped  silicon  has  yet  to  be  fully  characterized  at  high 
frequencies,  though  some  preliminary  results  have  been 
reported  [1].  Also,  metals  such  as  copper  and  aluminum 
have  long  demonstrated  high-frequency  conductivi¬ 
ties  significantly  lower  than  predicted  by  the  Drude 
model  [2]. 

The  Poisson  solver  incorporated  into  the  typical  En¬ 
semble  Monte  Carlo  (EMC)  implementation  assumes 
quasi- static  fields.  As  the  stimulating  frequency  ap¬ 
proaches  the  THz  regime,  this  assumption  loses  validity. 
Plasma  frequencies  and  characteristic  scattering  rates  of 
conductive  media  often  fall  within  the  THz  frequency 
range.  At  high  frequencies  fully  electrodynamic  field 
calculations  are  required  to  ensure  accuracy.  The  finite- 
difference  time-domain  (FDTD)  method  is  a  fully  elec¬ 
trodynamic,  highly  accurate  solver  based  on  Maxwell’s 
curl  equations  [3].  In  the  combined  EMC-FDTD  solver 
electric  and  magnetic  fields  from  the  FDTD  solver  influ¬ 
ence  EMC  carrier  motion.  Microscopic  currents  result¬ 
ing  from  carrier  motion  in  the  EMC  correspondingly  in¬ 
fluence  new  FDTD-computed  field  values.  A  simulation 
flowchart  is  shown  in  Fig.  1. 

Past  research  explored  the  combined  EMC-FDTD 
solver  in  the  context  of  AC  device  characterization  [4]. 
The  necessary  consideration  of  a  DC  component  requires 
continued  Poisson’s  equation  calculations  in  addition 
to  FDTD,  necessitating  substantial  computational  re¬ 
sources  and  limiting  solver  applicability. 

In  this  paper  we  examine  the  characteristics  of  the 
coupled  EMC-FDTD  solver  in  the  context  of  high- 
frequency  electromagnetic  wave  interactions  with  bulk 
doped  silicon.  The  Poisson  solver  is  not  necessary 
for  this  calculation.  We  characterize  the  EMC-FDTD 
solver  in  terms  of  grid  spacing,  carrier  ensemble  size, 
and  averaging  technique. 

We  use  two-dimensional  (2D)  computational  do¬ 
mains  for  both  FDTD  and  EMC.  The  FDTD  domain 
boundaries  are  treated  with  perfectly-matched  layer 
absorbing  boundary  conditions  for  conductive  media, 
allowing  finite-grid  represention  of  an  infinite  space. 
The  grid  is  defined  in  the  x  —  y  plane.  We  assume  a  TEZ 
mode  for  the  electromagnetic  wave;  thus  the  FDTD 
grid  comprises  Ex ,  Ey,  and  Hz  field  components.  The 
simulation  testbed  is  a  semi-infinite  half  space  of  doped 
silicon  (Regions  B  and  C  in  Fig.  2)  and  a  semi-infinite 
half  space  of  air  (Region  A  in  Fig.  2).  Region  A  is  as¬ 
signed  a  dielectric  constant  of  1  (and  zero  conductivity). 
A  dielectric  constant  of  11.8  is  assigned  to  Regions  B 


and  C.  Region  C  is  filled  with  an  assumed  value  for  the 
continuous  bulk  conductivity  of  doped  silicon,  (Jbuik- 
Region  B  is  in  continuous  interaction  with  the  EMC; 
the  conductivity  terms  in  the  FDTD  update  equations 
for  Ex  and  Ey  are  removed  from  Region  B,  and  interac¬ 
tion  of  FDTD  fields  and  EMC  currents  is  enforced.  We 
use  EMC  reflecting  boundary  conditions  at  the  vertical 
edges  of  Region  B  and  periodic  boundary  conditions  at 
the  horizontal  edges.  An  incident  0.5  THz  plane  wave  is 
introduced  via  the  analytic  field  propagation  total-field 
scattered-field  (AFP-TFSF)  formulation  [5]. 

Fig.  2  shows  the  2D  spatial  distribution  of  the  Ey 
phasor  amplitude,  extracted  over  several  periods  of  elec¬ 
tromagnetic  wave  oscillation.  The  field  amplitude  is 
attenuated  as  a  function  of  depth  in  both  Regions  B 
and  C;  this  macroscopic  phenomenon  is  indicative  of  the 
skin  effect  in  the  coupled  region.  The  effective  linear- 
regime  conductivity  observed  in  Region  B  is  calculated 
as  dobs  =  M{E  •  J*}/|E|2.  Significant  noise  in  \Ey\, 
resulting  from  electron  motion,  is  visible  in  Region  B. 
Noise  levels  in  Jx,  Jy,  and  Ex  in  that  region  are  simi¬ 
larly  impacted  by  electron  motion.  Accurate  calculation 
of  dobs  requires  careful  consideration  of  the  effects  of  this 
noise. 

Fig.  3  shows  the  observed  noise  of  \Ey\  in  Region 
B  for  several  configurations  of  grid  cell  size  Ax  and 
electron  ensemble  size  Ne.  Larger  carrier  ensembles 
show  decreased  noise,  as  expected.  This  improvement  is 
dwarfed,  however,  by  the  improvements  gained  through 
increased  grid  cell  size.  In  direct  contrast  to  the  EMC- 
Poisson  solver,  in  which  smaller  grid  cells  improve  ac¬ 
curacy,  the  EMC-FDTD  solver  produces  cleaner  results 
for  larger  grid  cells.  Instead  of  choosing  several  grid 
points  per  Debye  length,  as  per  Poisson  solver  accuracy 
and  stability  requirements,  the  EMC-FDTD  grid  cell 
size  may  be  defined  by  the  much  less  restrictive  FDTD 
numerical  dispersion  requirements,  which  only  require 
Ax  to  be  smaller  than  roughly  1  / 20th  of  the  smallest 
wavelength  of  interest. 

Noise  in  the  current  density  and  field  phasor s  may 
be  reduced,  and  more  accurate  conductivity  values  cal¬ 
culated,  through  spatial  averaging  of  extracted  phasor 
quantities  prior  to  the  aQ bs  calculation.  Fig.  4  shows  aQ bs 
plotted  with  respect  to  the  size  of  the  square  region  over 
which  the  phasor  quantities  are  averaged.  Noise  is  suffi¬ 
ciently  reduced  for  large  averaging  region  size  to  ensure 
convergence  in  cr0bs-  Fig-  4  shows  that  larger  ensembles 
also  allow  convergence  in  (Jobs,  supporting  the  assertion 
that  reduced  phasor  noise  permits  convergence  (Jobs- 

Finally,  in  Fig.  4  we  examine  the  importance  of 
impedance  matching  between  the  EMC-FDTD  coupled 


Region  B  and  the  surrounding  pure-FDTD  Region  C 
(filled  with  continuous  conductivity).  These  results 
show  dobs  to  be  insensitive  to  impedance  mismatch 
between  Regions  B  and  C. 

This  preliminary  work  on  a  combined  EMC-FDTD 
solver  demonstrates  the  feasibility  of  a  first-principles 
approach  to  predicting  the  effective  conductivity  of 
doped  semiconductors  at  THz  and  sub-THz  frequen¬ 
cies.  Further  work  is  needed  to  fully  evaluate  the 
performance  of  this  approach. 
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Figure  1:  Simulation  flowchart.  The  EMC  and  FDTD 
simulations  are  initialized  independently.  The  time  step 
counter  is  n.  Currents  and  fields  pass  between  the  two 
simulation  tools  at  every  time  step. 


Figure  2:  Amplitude  decay  of  Ey  field  phasor  extracted 
over  several  periods  of  electromagnetic  wave  oscillation. 
Air  fills  Region  A,  while  Regions  B  and  C  are  doped  sil¬ 
icon.  The  material  interface  is  indicated  by  the  vertical 
solid  white  line.  An  incident  plane  wave  is  sourced  from 
the  TFSF  boundary,  shown  by  the  dashed  line.  Region 
B  is  the  EMC-FDTD  coupled  region,  whose  boundary 
is  marked  by  a  solid  white  box,  and  Region  C  is  pure- 
FDTD  doped  silicon. 


Figure  3:  Average  noise  of  Ey  in  the  EMC-FDTD  cou¬ 
pled  region  for  several  values  of  electron  ensemble  size 
and  grid  cell  size,  where  Axref  ~  400  nm  is  a  reference 
grid  cell  size.  Ey°lse  is  calculated  as  {\Ey\2}  —  \(Ey)\2 
and  is  normalized  by  \Ey  \  at  the  material  interface.  In¬ 
creasing  the  size  of  the  ensemble  reduces  noise,  as  ex¬ 
pected,  but  increasing  the  grid  cell  size  produces  a  much 
larger  improvement. 


Figure  4:  Observed  conductivity  with  varied  averaging 
region  size,  electron  ensemble  size,  and  surrounding  bulk 
conductivity.  No  is  a  reference  ensemble  size,  typically 
O(105).  Both  electron  ensemble  size  and  averaging  re¬ 
gion  size  increase  cr0bs-  ^obs  is  insensitive  to  impedance 
mismatch  between  Regions  B  and  C. 

[1]  T.  I.  Jeon  and  D.  Grischkowsky,  Phys.  Rev.  Lett. 
Vol.  78,  pp.  1106-1109,  1997. 

[2]  A.  B.  Pippard,  Proc.  R.  Soc.  A-Math.  Phys.  Eng. 
Sci.  Vol.  191,  pp.  385-399,  1947. 

[3]  A.  Taflove  and  S.  C.  Hagness,  Computational  Elec¬ 
trodynamics:  The  Finite- Difference  Time-Domain 
Method ,  Artech  House,  2nd  ed.,  2000. 

[4]  J.  S.  Ayubi-Moak,  S.  M.  Goodnick,  and  M.  Saran- 
iti,  J.  Comput.  Electron.  Vol.  5,  pp.  415-418,  2006. 

[5]  J.  B.  Schneider,  IEEE  Trans.  Antennas  Propag. 
Vol.  52,  pp.  3280-3287,  2004. 


Dielectric  Characterization  of  Carbon  Nanotube  Contrast  Agents  for 
Microwave  Breast  Cancer  Detection 


Alireza  Mashal*(I),  Balaji  Sitharaman  (2),  John  H.  Booske(I),  and  Susan  C.  Hagness(I) 

(1)  Department  of  Electrical  Engineering,  University  of  Wisconsin,  Madison,  WI  53706 
(2)  Department  of  Biomedical  Engineering,  Stony  Brook  University, 

Stony  Brook,  NY  11794 

E-mail:  amashal@wisc.edu,  booske@engr.wisc.edu,  hagness@engr.wisc.edu 

Abstract 

We  are  investigating  the  feasibility  of  single-walled  carbon  nanotubes  (SWCNTs)  as 
contrast  agents  for  microwave  breast  cancer  detection.  We  hypothesize  that  the 
accumulation  of  SWCNTs  in  a  tumor  will  enhance  the  dielectric  contrast  between  normal 
and  malignant  tissue  and  potentially  improve  the  efficacy  of  microwave  imaging 
techniques.  As  an  initial  step,  we  constructed  tissue-mimicking  phantom  materials  with 
varying  concentrations  of  SWCNTs.  We  characterized  their  dielectric  properties  and 
estimated  the  pressure  wave  induced  in  the  phantom  materials  under  microwave 
illumination.  At  SWCNT  concentrations  of  less  than  0.5%  by  weight,  we  observed 
significant  changes  in  both  the  relative  permittivity  and  effective  conductivity  of  these 
SWCNT/tissue-mimicking  material  mixtures.  Our  estimates  suggest  a  similarly 
significant  change  in  the  thermoacoustic  response  relative  to  the  phantom  material 
mixtures  without  carbon  nanotubes. 


1.  Introduction 

Microwave  detection  of  breast  cancer  has  been  under  intense  investigation  for  the  past 
decade  because  of  its  potential  as  a  low  cost,  non-ionizing  modality  that  offers  3D 
tomographic  imaging  capabilities.  Until  recently,  work  in  this  area  has  primarily  focused 
on  exploiting  the  intrinsic  dielectric  contrast  between  normal  and  malignant  tissue. 
However,  a  published  large-scale  dielectric  spectroscopy  study  conducted  by  the 
Universities  of  Wisconsin-Madison  and  Calgary  has  shown  that  while  the  contrast  in  the 
microwave-frequency  dielectric  properties  between  malignant  tissues  and  healthy 
adipose-dominated  tissues  in  the  breast  is  as  large  as  10:1,  the  contrast  between  malignant 
and  normal  glandular/fibroconnective  tissues  in  the  breast  is  no  more  than  about  10%  [1]. 
Since  most  tumors  originate  in  glandular  tissue,  these  results  suggest  that  solely 
exploiting  the  intrinsic  dielectric  contrast  between  normal  and  malignant  tissue  will  result 
in  a  more  challenging  cancer  detection  scenario  than  previously  thought. 

We  have  recently  initiated  an  investigation  of  exogenous  contrast  agents  to  determine 
how  they  affect  the  sensitivity  of  microwave  imaging  modalities.  In  [2,  3],  we  describe 
the  dielectric  contrast-enhancement  effects  of  gas-filled  microbubbles  and  metallic 
nanoparticles,  while  in  [4]  we  describe  the  impact  of  microbubbles  on  3D  microwave 
tomographic  breast  imaging.  Here  we  propose  the  use  of  single  wall  carbon  nanotubes 
(SWCNTs)  as  a  contrast  agent  for  microwave  imaging.  Biocompatible  carbon  nanotubes 
are  being  broadly  explored  for  a  variety  of  diagnostic  and  therapeutic  applications  in 
medicine  such  as  nanoscale  biosensors,  drug  delivery  for  cancer  treatment,  novel 
biomaterials,  and  molecular  imaging  [5].  SWCNTs  that  accumulate  in  tumors  via 
passive  or  active  mechanisms  could  alter  the  effective  complex  permittivity  of  the  tumor 
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when  compared  to  a  tumor  in  its  pre-contrast  stage.  This  change  in  dielectric  properties 
can  be  exploited  via  differential  imaging  to  improve  sensitivity.  In  addition,  the  increase 
in  effective  conductivity  could  be  advantageous  in  thermal  therapy  applications.  To  test 
the  feasibility  of  these  carbon  nanotubes  as  a  contrast  agent  for  microwave  applications, 
we  characterize  the  dielectric  properties  of  tissue  mimicking  (TM)  materials  mixed  with 
various  concentrations  of  SWCNTs.  We  then  evaluate  the  impact  of  the  dielectric 
contrast  enhancement  on  microwave-induced  thermoacoustic  imaging  by  estimating  the 
amplitude  of  the  induced  pressure  wave  in  a  thermoacoustic  imaging  system. 

2.  Materials  and  methods 

Lazebnik  et  al  [6]  discusses  the  construction  of  TM  materials  that  approximate  the 
dielectric  properties  of  a  variety  of  human  soft  tissues  over  the  microwave  frequency 
range  that  spans  0.5  GHz  to  20  GHz.  The  results  reported  in  [1]  show  that  a  10%  oil- 
gelatin  mixture  adequately  describes  the  microwave  properties  of  malignant  breast  tissue. 
Therefore,  for  this  study  we  constructed  10%  oil-gelatin  TM  materials  by  following  the 
procedure  outlined  in  [6]  and  adding  the  step  of  introducing  of  SWCNTs.  We  mixed  in 
the  SWCNTs  by  first  adding  the  indicated  amount  of  SWCNTs  to  a  1%  by  weight 
mixture  of  Pluronic  (FI 27)  deionized  water  mixture.  The  SWCNTs  we  used  had  an  outer 
diameter  of  1-2  nm,  length  of  5-30  |Lim,  and  purity  of  greater  than  90%  by  weight.  This 
mixture  was  then  sonicated  using  a  probe  sonicator  for  twenty  minutes  at  a  120  W  power 
level.  We  made  various  concentrations  (1  mg/mL,  2  mg/mL,  and  3  mg/mL)  of  SWCNT 
in  deionized  water  and  substituted  these  mixtures  for  the  pure  deionized  water  in  the  TM 
recipe  described  in  [6].  We  poured  the  liquid  oil-gelatin  mixtures  into  small  covered 
glass  jars  and  allowed  the  mixture  to  gel  unperturbed. 

We  characterized  the  dielectric  properties  of  the  mixtures  using  a  wideband  open-ended 
coaxial  probe  technique  [7],  We  positioned  the  tip  of  the  probe  against  the  surface  of  the 
TM  mixture,  being  careful  not  to  compress  the  gel  excessively  or  puncture  the  surface  as 
described  in  [6].  We  made  measurements  on  three  different  samples  of  each  TM  mixture 
and  averaged  the  results. 


3.  Results  and  discussion 

Figures  1(a)  and  (b)  show  the  measured  relative  permittivity  and  effective  conductivity  of 
the  pure  TM  material  and  the  TM  materials  with  three  different  SWCNT  concentrations 
over  a  frequency  range  of  0.6  -  20  GHz.  Each  wideband  curve  represents  the  average  of 
measurements  on  three  samples  of  the  given  mixture  and  the  vertical  bars  depict  the 
deviation  of  the  individual  measurements  from  the  average  value.  Table  1  summarizes 
this  data  at  3  GHz  and  quantifies  the  percent  change  in  dielectric  properties  due  to  the 
introduction  of  the  SWCNTs.  Our  results  show  that  small  increases  in  the  concentration 
of  SWCNTs  lead  to  a  significant  increase  in  both  the  relative  permittivity  and  effective 
conductivity  of  the  TM  material. 


(a)  (b) 

Figure  1.  (a)  Relative  permittivity  and  (b)  effective  conductivity  of  TM  mixture  with 
varying  concentrations  of  SWCNTs. 


Table  L  Summary  of  dielectric  spectroscopy  data  at  3  GHz. 


Concentration  of 
SWCNT  in  water 

%  by  weight  of 
SWCNTs 

Average  f, 

(3  GHz) 

Average  a 
(3  GHz) 

Percent 
change  % 

Percent 
change  a 

0 

0 

53.22 

2.31  S/m 

- 

- 

1  mg/mL 

0.074  % 

57.99 

2.69  S/m 

9% 

16% 

2  mg/mL 

0.148% 

63.59 

3.19  S/m 

19% 

38% 

3  mg/mL 

0.222  % 

72.75 

4.18  S/m 

37% 

81  % 

The  changes  in  dielectric  properties  can  lead  to  an  improvement  in  the  sensitivity  of 
microwave  imaging  techniques,  such  as  microwave-induced  thermoacoustic  imaging 
(MI-TAT).  MI-TAT  is  a  hybrid  imaging  modality  that  exploits  dielectric-properties 
contrasts  at  microwave  frequencies  while  creating  images  with  ultrasound-quality 
resolution.  In  MI-TAT,  the  tissue  is  irradiated  with  a  sub-microsecond  electromagnetic 
pulse,  which  is  then  selectively  absorbed  by  the  higher-conductivity  tissues  and  induces 
acoustic  waves.  These  acoustic  waves  are  detected  by  ultrasound  transducers  and 
processed  for  image  reconstruction.  We  estimate  the  thermoacoustic  response  of  a 
simple  target  by  calculating  the  amplitude  of  the  induced  pressure  wave  assuming  stress 
confinement,  uniform  heating,  and  neglecting  thermal  diffusion  [8]: 


AP  =  B/3e-AT  =B-J3e 


SAR  t_B_/3^t  o\E\2 


2  P 


oca. 


(i) 


Here,  /?e  is  the  average  volumetric  expansion  coefficient  over  small  temperature  ranges,  B 
is  the  bulk  modulus,  A T  is  the  change  in  local  temperature,  r  is  the  microwave  pulse 
length,  Cp  is  the  specific  heat.  SAR  is  the  specific  absorption  rate,  defined  as  o\Ef/2p, 
where  E  is  the  electric  field,  a  is  the  effective  electrical  conductivity,  and  p  is  the  mass 
density. 

Equation  1  shows  that  the  amplitude  of  the  induced  pressure  wave  is  proportional  to  the 
effective  conductivity  of  the  tissue.  If  we  assume  that  the  carbon  nanotubes  do  not 
significantly  alter  the  thermomechanical  properties  of  the  tissue  or  the  electric  field 
intensity  penetrating  the  tumor,  then  we  can  infer  that  the  dielectric  contrast  provided  by 
the  SWCNT’s  correlate  to  a  similar  contrast  in  the  thermoacoustic  response  and  hence 
can  aid  the  sensitivity  of  MI-TAT  imaging. 


Conclusions 


We  characterized  the  dielectric  properties  of  TM  materials  with  varying  concentrations  of 
SWCNTs.  Our  results  indicate  that  concentrations  as  small  as  0.222%  by  weight  of 
SWCNTs  alter  the  relative  permittivity  of  the  TM  phantom  material  by  37%  and  the 
effective  conductivity  by  81%  when  compared  to  the  TM  material  without  SWCNTs. 
Our  estimates  show  a  similar  change  in  the  induced  pressure  in  a  thermoacoustic  imaging 
system.  These  changes  in  dielectric  properties  and  thermoacoustic  response  at  such  low 
concentrations  of  carbon  nanotubes  are  promising  and  suggest  that  SWCNTs  have  the 
potential  to  be  effective  as  a  contrast  agent  in  microwave  imaging  modalities. 
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